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ABSTRACT 

Significant  advances  have  been  made  in  several  aspects  of  the  computational  modelling  of 
turbulent  combustion.  PDF  model  calculations  have  been  performed  of  turbulent  piloted-jet  non- 
premixed  flames.  The  results  demonstrated  the  ability  of  the  methodology  to  account,  accurately, 
for  the  local  extinction  and  reignition  observed  experimentally  in  these  flames.  It  was  shown 
that  these  flames  can  be  sensitive  to  the  temperature  of  the  pilot  and  to  radiative  heat  loss.  A 
new  approach  has  been  developed  for  the  efficient  computational  implementation  of  combustion 
chemistry.  The  rate-controlled  constrained  equilibrium  method  has  been  combined  with  the  in 
situ  adaptive  tabulation  algorithm  to  produce  a  unified  dimension-reduction/storage-retrieval 
methodology  for  the  computationally-efficient  implementation  of  combustion  chemistry.  Test 
calculations  demonstrated  that  this  methodology  has  comparable  accuracy  to  augmented  reduced 
mechanisms.  Ideas  from  the  conditional  moment  closure  and  the  mapping  closure  have  been 
combined  to  produce  a  new  approach  for  modeling  molecular  mixing  in  turbulent  reactive  flows. 
The  new  methodology  has  been  shown  to  describe  accurately  (for  the  first  time)  the  mixing  of 
two  scalars.  A  methodology  has  been  developed  for  obtaining  stochastic  models  for  Lagrangian 
velocity  and  acceleration  based  on  DNS  data  from  homogeneous  turbulent  shear  flow.  It  has 
been  shown  that  the  acceleration  model  provides  a  remarkably  accurate  representation  of  the 
observed  Lagrangian  velocity-acceleration  two-time  correlations.  In  collaboration  with  the  group 
of  Prof.  P.  Givi,  advances  have  been  made  in  the  implementation  of  a  combined  LES/PDF 
methodology  for  modeling  turbulent  reactive  flows.  The  approach  based  on  the  velocity  filtered 
density  function  has  been  applied  to  a  spatially-developing  mixing  layer  and  shown  to  account 
well  for  the  major  processes  in  this  flow. 

INTRODUCTION 

The  design  process  for  gas-turbine  combustors  and  aerospace  propulsion  systems  could  be 
significantly  improved  if  accurate  and  affordable  CRD  tools  were  available.  While  turbulent 
combustion  models  are  used  in  the  design  process,  the  models  currently  employed  are  not 
sufficiently  accurate.  PDF  methods  promise  the  capability  of  greater  accuracy,  through  their 
ability  to  treat  the  chemistry  in  sufficient  detail,  and  to  fully  account  for  turbulence-chemistry 
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interactions.  The  work  performed  in  this  research  project  has  significantly  contributed  to  the 
development  and  demonstration  of  different  aspects  of  PDF  methods. 

The  research  has  focused  on  the  following  topics. 

1/  PDF  calculations  of  turbulent  non-premixed  piloted  jet  flames. 

2/  The  development  of  a  new  methodology  (based  on  rate-controlled  constrained  equilibrium 
and  in  situ  adaptive  tabulation)  for  the  efficient  implementation  of  combustion  chemistry. 

3/  The  development  of  a  new  approach  (related  to  the  mapping  closure  and  the  conditional 
moment  closure)  for  modeling  mixing  in  turbulent  reactive  flows. 

4/  The  development  of  stochastic  Lagrangian  models  for  velocity  and  acceleration  in  turbulent 
flows,  and  a  methodology  to  determine  the  model  coefficients  involved  from  DNS  data. 

5/  The  implementation  and  demonstration  of  the  approach  based  on  the  combination  of  large 
eddy  simulation  and  PDF  methods. 

The  work  on  each  topic  is  described  in  the  following  sections,  and  more  completely  in  the 
publications  given  below. 

PDF  CALCULATIONS  OF  PILOTED  NONPREMIXED  FLAMES 

The  piloted  nonpremixed  flames  studied  experimentally  at  Sandia  (Barlow  &  Frank  1998) 
provide  an  excellent  test  of  turbulent  combustion  models.  These  flames  show  distinct  levels  of 
interaction  between  turbulence  and  chemistry  because  of  the  increasing  jet  bulk  velocities  from 
flame  D  to  F:  flame  D  is  close  to  equilibrium  with  a  small  amount  of  local  extinction,  whereas 
flame  F  is  on  the  verge  of  global  extinction.  In  each  of  these  flames,  the  amount  of  local 
extinction  reaches  a  peak  at  an  axial  distance  of  about  30  jet  radii,  with  re-ignition  occurring 
downstream.  Several  advanced  approaches  based  on  LES,  CMC  and  PDF  methods  have  been 
applied  to  compute  these  phenomena  and  have  made  significant  progress.  Notably,  the  joint 
PDF  calculations  of  these  flames  by  Xu  and  Pope  (2000)  and  Lindstedt  et  al.  (2000)  show  the 
best  detailed  agreement  obtained  between  computations  and  the  experimental  data. 

The  PDF  calculations  of  Xu  &  Pope  (2000)  and  Tang  et  al.  (2000)  are  capable  of  calculating, 
quantitatively,  the  observed  phenomena  of  local  extinction  and  reignition.  These  calculations  are 
based  on  the  modelled  transport  equation  for  the  joint  PDF  of  velocity,  turbulence  frequency,  and 
composition.  The  sub-models  of  this  method  include  the  simplified  Langevin  model  (SLM)  for 
velocity  and  the  Jayesh-Pope  model  (JPM)  for  turbulent  frequency  (see  Pope  2000).  The 
molecular  mixing  is  modeled  by  the  Euclidean  minimal  spanning  tree  (EMST)  model  of 
Subramaniam  &  Pope  (1998),  which  features  mixing  locally  in  the  composition  space  through 
interacting  particles  with  neighboring  particles.  The  reaction  mechanism  used  is  the  19  species, 
15-step  augmented  reduced  mechanism  of  Sung  et  al.  (1998)  which  includes  NO  chemistry,  and 
is  denoted  by  ARM2.  The  chemical  reaction  calculations  are  performed  using  the  in  situ  adaptive 
tabulation  (IS AT)  algorithm  (Pope  1997).  It  should  be  pointed  out  that  for  each  full-scale  PDF 
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method  calculation,  the  solution  to  the  reaction  equation  system  (20  dimensional)  is  required  (of 
order)  109  times.  IS  AT  can  handle  these  computations  economically  and  accurately.  Recent 
work  on  PDF  methods  for  these  flames-now  described-concerns  sensitivity  to  the  pilot  flame 
temperature  and  to  radiative  heat  loss. 

It  was  observed  by  Xu  &  Pope  (2000)  that  calculations  of  flame  F  exhibit  some  sensitivity  to  the 
pilot  temperature  Tp  which  is  specified  as  a  boundary  condition.  The  experimental  data  show  Tp 
in  the  range  1860K-1880K,  but  the  experimental  accuracy  may  be  no  greater  than  10-20K.  This 
influence  was  studied  systematically  by  performing  calculations  of  flames  D  and  F  with  pilot 
temperatures  of  Tp  =  1860,  1870  and  1880K.  For  flame  D  it  is  found  that  the  calculations  are 
insensitive  to  Tp.  But,  as  shown  on  Fig.  1,  flame  F  exhibits  extreme  sensitivity.  For  example,  at 
x/Rj  =  15  and  30  the  peak  temperature  decreases  by  about  500K  and  the  mass  fractions  of  OH 
decrease  by  a  factor  of  more  than  two  with  a  10K  decrease  in  pilot  temperature.  Similar  trends 
exist  for  other  variables.  Particularly  for  NO,  the  results  calculated  using  lower  pilot 
temperatures  give  a  perfect  match  with  the  experimental  data  at  the  first  two  locations  shown  in 
the  figure.  However,  further  downstream,  all  the  modeled  NO  profiles  overshot  the  peak  value  by 
a  factor  of  more  than  two  and  do  not  predict  correctly  on  the  fuel  rich  side,  although  the 
temperature  profiles  seem  to  be  satisfactory. 

Figure  2  shows  another  manifestation  of  the  sensitivity  to  the  pilot  temperature.  The  burning 
index  (BI)  is  defined  to  be  unity  for  complete  combustion  and  zero  for  complete  extinction.  The 
figure  shows  substantially  different  results  for  Tp  =  1860K  and  Tp  =  1880K,  with  the 
experimental  data  generally  falling  between  these  two  calculated  values. 

The  effects  of  radiative  heat  loss  were  investigated  by  performing  "adiabatic"  and  "radiant" 
calculations.  In  the  former  all  heat  loss  is  neglected:  in  the  latter,  radiative  heat  loss  is  accounted 
for  from  the  primary  radiating  species,  C02,  H20,  CO  and  CH4.  We  adopt  an  optically-thin  limit 
radiation  model,  although  the  validity  of  this  model  for  the  4.3-micron  band  of  C02  is  still  in 
debate.  Implemented  in  the  framework  of  ISAT,  the  model  includes  the  above  four  gas-phase 
emitting  species  with  their  Plank  mean  absorption  coefficients  calculated  by  RADCAL. 

For  flame  D  (not  shown)  the  effect  of  radiation  is  to  reduce  the  peak  temperature  by  about  30K 
and  to  decrease  the  peak  NO  by  15%.  Figure  3  presents  the  conditional  mean  profiles  of  four 
scalars,  and  shows  completely  different  picture  from  the  flame  D  results.  For  Tp  =  1880K,  the 
inclusion  of  radiation  induces  significant  decreases  in  temperature  and  species  mass  fractions  at 
the  first  three  axial  locations.  The  largest  differences  appear  at  x/Rj  =  30  where  the  peak 
temperature  decreases  more  than  500K  and  the  species  mass  fractions  decrease  by  a  factor  of  two 
or  three.  This  fact  indicates  that  thermal  radiation  can  significantly  alter  the  local  extinction 
status  in  this  flame:  not  only  is  the  NO  chemistry  strongly  influenced  by  radiation,  but  also  the 
reactions  of  other  species  such  as  OH  and  CO.  The  last  column  of  Fig.  3  tells  us  that  further 
downstream,  the  flame  becomes  closer  to  the  equilibrium  state  as  re-ignition  takes  place  and  the 
radiation  tends  to  be  less  important. 

Evidently,  flame  F  is  extremely  sensitive  to  a  decrease  in  temperature,  whether  it  arises  from  Tp 
or  from  radiation. 
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EFFICIENT  IMPLEMENTATION  OF  COMBUSTION  CHEMISTRY 


It  is  computationally  prohibitive  to  use  detailed  hydrocarbon  chemistry  directly  in  turbulent 
combustion  calculation.  Two  separate  approaches  have  been  taken  to  reduce  the  computational 
burden:  dimension  reduction,  and  storage/retrieval.  Tang  &  Pope  (2002)  have  combined  these 
two  approaches  into  a  unified  methodology.  Dimension  reduction  is  achieved  through  rate- 
controlled  constrained  equilibrium  ((RCCE  Keck,  1990);  and  storage/retrieval  through  the  IS  AT 
algorithm  Pope  (1997)).  In  this  context,  RCCE  is  preferred  over  other  reduction  methodologies 
(e.g.,  QSSA,  ILDM),  because  of  the  guaranteed  existence  and  continuity  of  the  implied  low¬ 
dimensional  manifold. 

The  combined  ISAT-RCCE  methodology  is  tested  for  a  pairwise-mixing  stirred  reactor  (PMSR) 
using  the  31 -species  GRI  1.2  mechanism  for  methane.  Three  different  tests  (referred  to  as  Ci,  C2 
and  C3)  are  performed.  In  C\  the  constrained  species  are  H2O,  CO2,  O2,  CH4  and  CO.  Three  more 
species  (H2,  OH,  and  O )  are  added  to  form  the  constraint  subspace  in  C2,  and  in  C3  another  three 
species  ( C//3 ,  C2H2  and  C2H4)  are  included.  Additional  constraints  are  imposed  on  all  4  elements 
and  enthalpy,  and  hence  the  dimension  reductions  are  from  32  to  10,  13  and  16  for  the  three 
cases,  respectively.  To  test  the  accuracy  of  the  algorithm,  we  solve  the  entire  ODE  system  (32- 
dimensional)  by  direct  integration  (DI)  to  get  the  accurate  solution.  Figure  4  shows  the  relative 
error  in  species  compositions,  temperature  and  density  against  their  reference  values  for  one 
particle  (advanced  over  2000  time  steps  in  the  statistically  stationary  state).  It  can  be  observed 
that,  as  the  number  of  constraints  increases,  the  error  in  the  constrained  species  decreases.  For 
C3,  the  relative  errors  in  major  species  (including  CO  and  H2)  are  under  3%  with  the  errors  of 
other  constrained  species  being  less  than  10%. 

In  Figure  5,  the  accuracy  of  ISAT-RCCE  is  compared  to  that  of  the  augmented  reduced 
mechanism  of  Sung  et  al.  (1998)-which  is  based  on  the  same  detailed  mechanism  and  which  has 
the  same  number  of  degrees  of  freedom.  It  may  be  seen  that  the  two  methods  have  comparable 
accuracy. 

MODELLING  TURBULENT  MIXING 

In  PDF  methods  for  turbulent  combustion,  the  modeling  of  molecular  diffusion  is  both  crucial 
and  difficult.  In  Klimenko  &  Pope  (2002)  a  new  methodology  is  developed-multiple  mapping 
conditioning  (MMC)-which  combines  ideas  from  the  mapping  closure  (Chen  et  al.  1989),  and 
from  the  conditional  moment  closure  (Klimenko  &  Bilger  1999).  In  part,  this  approach  extends 
the  particle  implementation  of  the  mapping  closure  to  multiple  scalars. 

Remarkably,  the  MMC  model  admits  an  analytic  solution  for  the  case  of  two  passive  scalars 
evolving  from  a  triple-delta-function  initial  condition.  This  case  was  studied  using  DNS  by 
Juneja  &  Pope  (1996),  with  the  three  delta  functions  located  at  the  vertices  of  an  equilateral 
triangle  in  the  two-dimensional  composition  space.  The  evolution  predicted  by  MMC  (Fig.  6)  is 
in  excellent  agreement  with  the  DNS.  No  other  mixing  model  has  been  shown  to  be  even 
qualitatively  correct  for  this  case. 
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STOCHASTIC  MODELLING  OF  VELOCITY  AND  ACCELERATION 

In  PDF  methods,  the  turbulence  modeling  is  embodied  in  a  stochastic  model  for  the  velocity 
following  a  fluid  particle  (see  e.g.,  Pope  2000).  The  standard  model-the  generalized  Langevin 
model-involves  tensor  coefficients.  In  Pope  (2002a)  a  methodology  is  developed  to  determine 
these  coefficients  from  DNS  data.  In  Pope  (2002b)  this  methodology  is  extended  to  a  stochastic 
model  for  acceleration,  which  is  a  natural  way  to  incorporate  Reynolds-number  effects. 

Figure  7  shows  the  velocity-acceleration  autocovariances  predicted  by  the  model  compared  to 
the  DNS  data  of  Sawford  &  Yeung  (2000).  As  may  be  seen,  the  model  is  able  to  provide  an 
accurate  representation  of  these  fundamental  statistics. 

LARGE  EDDY  SIMULATION 

The  PDF  calculations  reported  above  (e.g.,  Xu  &  Pope  2000)  are  based  on  a  completely 
statistical  approach.  For  some  turbulent  reacting  flow,  especially  those  with  large-scale  unsteady 
motions,  there  is  good  reason  to  use  large  eddy  simulation  (LES).  In  this  approach,  the  large 
scales  are  treated  deterministically,  and  the  small  scales  statistically.  It  is  important  to  appreciate 
that  for  turbulent  combustion,  the  subgrid  scale  modeling  pertaining  to  reaction  required  in  LES 
is  similar  to,  and  just  as  crucial  as,  the  modeling  of  turbulence-combustion  interactions  in 
statistical  approaches.  It  is  natural,  therefore,  to  combine  the  ability  of  LES  to  represent  large- 
scale  unsteady  motions,  with  the  benefits  of  PDF  methods  for  modeling  the  turbulent-chemistry 
interactions.  The  PI  has  a  long-term  collaboration  with  the  group  of  Prof.  Peyman  Givi  to 
develop  this  combined  LES/PDF  methodology. 

Following  the  earlier  work  of  Colucci  et  al.  (1998)  and  Jaberi  et  al.  (1999),  recently  Gicquel  et 
al.  (2002)  have  extended  the  methodology  to  consider  the  velocity  filtered  density  function 
(VFDF).  In  this  methodology,  the  effects  of  the  unresolved  subgrid  scales  (SGS)  are  taken  into 
account  by  considering  the  joint  probability  density  function  of  all  of  the  components  of  the 
velocity  vector.  An  exact  transport  equation  is  derived  for  the  VFDF  in  which  the  effects  of  the 
SGS  convection  appear  in  closed  form.  The  unclosed  terms  in  this  transport  equation  are 
modeled.  A  system  of  stochastic  differential  equations  (SDEs)  which  yields  statistically 
equivalent  results  to  the  modeled  VFDF  transport  equation  is  constructed.  These  SDEs  are 
solved  numerically  by  a  Lagrangian  Monte  Carlo  procedure  in  which  the  Ito  character  of  the 
SDEs  is  preserved.  The  consistency  of  the  proposed  SDEs  and  the  convergence  of  the  Monte 
Carlo  solution  are  assessed  by  comparison  with  results  obtained  by  an  Eulerian  LES  procedure  in 
which  the  corresponding  transport  equations  for  the  first  two  SGS  moments  are  solved.  The 
VFDF  results  are  compared  with  those  obtained  via  several  existing  SGS  closures.  These  results 
are  also  analyzed  via  a  priori  and  a  posteriori  comparisons  with  results  obtained  by  direct 
numerical  simulation  of  an  incompressible,  three-dimensional,  temporally  developing  mixing 
layer. 
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FIGURES 


Fig.  1.  Radial  profiles  of  mean  temperature  and  mass  fractions  in  flame  F:  sensitivity  to  pilot  temperature, 
Tp  Symbols,  experimental  data;  solid  line  Tp  —  1880K;  dashed  line  Tp  =  1870K;  dot-dashed  line  Tp  — 
1860K. 
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Fig.  3.  Effect  of  radiative  heat  loss  on  means  conditional  on  mixture  fraction  in  flame  F.  Symbols 
experimental  data;  solid  line,  adiabatic  calculation;  dashed  line,  radiant  calculation. 


Fig.4:  Relative  errors  for  different  species,  temperature  and  density,  against  the  reference  value  for  each 
quantity,  for  different  numbers  of  represented  species  in  ISAT-RCCE.  From  the  PMSR  test  case  of  Tang  & 
Pope  (2002). 
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Fig.  5:  Relative  errors  for  different  species,  temperature  and  density,  against  the  reference  value  for  each 
quantity,  from  the  PMSR  test  case:  comparison  of  ISAT-RCCE  (Tang  &  Pope  2002)  with  the  augmented 
reduced  mechanism  (Sung  et  al.  1998). 
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Fig.  6:  Temporal  evolution  (a-f)  of  the  joint  PDF  of  two  passive  scalars  from  a  triple-delta-function  initial 
condition  according  to  the  MMC  closure  of  Klimenko  &  Pope  (2002). 


Fig.  7:  Velocity-acceleration  Lagrangian  covariances:  symbols,  from  the  DNS  data  of  Sawford  &  Yeung 
(2001);  lines,  from  the  stochastic  model  of  Pope  (2002b). 
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Stochastic  Lagrangian  models  of  velocity  in  homogeneous 
turbulent  shear  flow 
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Stochastic  Lagrangian  models  for  the  velocity  following  a  fluid  particle  are  used  both  in  studies  of 
turbulent  dispersion  and  in  probability  density  function  (PDF)  modeling  of  turbulent  flows.  A 
general  linear  model  is  examined  for  the  important  case  of  homogeneous  turbulent  shear  flow,  for 
which  there  are  recent  direct  numerical  simulation  (DNS)  data  on  Lagrangian  statistics.  The  model 
is  defined  by  a  drift  coefficient  tensor  and  a  diffusion  tensor,  and  it  is  shown  that  these  are  uniquely 
determined  by  the  normalized  Reynolds-stress  and  timescale  tensors  determined  from  DNS.  With 
the  coefficients  thus  determined,  the  model  yields  autocorrelation  functions  in  good  agreement  with 
the  DNS  data.  It  is  found  that  the  diffusion  tensor  is  significantly  anisotropic — contrary  to  the 
Kolmogorov  hypotheses  and  conventional  modeling — which  may  be  a  low-Reynolds-number  effect. 
The  performance  of  two  PDF  models  is  also  compared  to  the  DNS  data.  These  are  the  simplified 
Lagrangian  model  and  the  Lagrangian  isotropization  of  production  model.  There  are  significant 
differences  between  the  autocorrelation  functions  generated  by  these  models  and  the  DNS  data. 
©  2002  American  Institute  of  Physics.  [DOI:  10.1063/1.1465421] 


I.  INTRODUCTION 

Homogeneous  turbulent  shear  flow  is  of  fundamental 
importance  in  the  development  of  models  for  inhomoge¬ 
neous  turbulent  flows.  Both  experiments1  and  direct  numeri¬ 
cal  simulations  (DNS)2  of  homogeneous  shear  flow  have 
been  performed  in  which  Eulerian  statistics  of  the  turbulence 
have  been  measured.  More  recently,  a  series  of  DNS  studies 
has  been  performed3-5  in  which  Lagrangian  statistics  have 
been  obtained  by  tracking  a  large  number  of  fluid  particles. 
These  studies  clearly  have  direct  relevance  to  stochastic  La¬ 
grangian  models6  of  turbulence,  which  model  the  motion  of 
fluid  particles  as  diffusion  processes  (i.e.,  continuous  Markov 
processes).7  The  purpose  of  this  paper  is  to  show  the  connec¬ 
tion  between  the  Lagrangian  velocity  autocovariance  tensor 
obtained  from  DNS  and  stochastic  Lagrangian  models  for 
fluid  particle  velocity. 

Stochastic  Lagrangian  models  for  the  velocity  of  a  fluid 
particle  arise  in  two  different  contexts:  turbulent 
dispersion;8-10  and  probability  density  function  (PDF) 
models. 11,12,7  In  both  cases  the  general  form  of  the  models 
considered  (when  applied  to  homogeneous  turbulence)  can 
be  written  as  the  linear  stochastic  differential  equation  (SDE) 

du^-AijUjdt+BijdWj,  (1) 

where  du(t)^u(t+dt)  —  u(t)  is  the  infinitesimal  increment 
of  the  fluctuating  component  of  velocity  u(r)  following  the 
fluid  particle;  we  refer  to  A (t)  as  the  drift  tensor;  B(r)  is  the 
diffusion  coefficient;  and  dW(t)  is  the  infinitesimal  incre¬ 
ment  of  a  vector- valued  Wiener  process  which  has  the  prop¬ 
erties  (dW)  =  0,  ( dW i  dWj)=dtSij .  Different  models  corre¬ 
spond  to  different  specifications  of  the  drift  tensor  A(t)  and 
diffusion  coefficient  B(r). 


For  statistically  stationary,  homogeneous  isotropic  turbu¬ 
lence  (with  no  mean  velocity  gradients)  the  only  sensible 
choice  of  coefficients  is 

Su 

A ij  yl  ’  ® 

and 


2m' 


/  2\  1/2 


(3) 


where  Th  is  the  Lagrangian  integral  timescale  and  uf  is  the 
turbulence  intensity  (i.e.,  the  rms  velocity  fluctuation).  Then, 
Eq.  (1)  reduces  to  an  independent  Lange vin  equation  for 
each  component  of  velocity 

dt  I  2u,2\  1/2 

du^-Ui—  +  |  —  I  dWt.  (4) 

This  model  dates  back  to  Taylor’s  1921  original  paper  on 
turbulent  dispersion.8  The  autocorrelation  function  given  by 
Eq.  (4)  is 

p(i)=(«1(f)«I(/+i))/M'2=exp(-|j|/rt),  (5) 

which  agrees  well  with  DNS  data13  (except  at  small  values  of 

M/rL). 

The  central  issue  addressed  here  is  the  appropriate  speci¬ 
fication  of  A  and  B  in  homogeneous  turbulent  shear  flow. 
This  has  been  considered  in  the  context  of  turbulent  disper¬ 
sion  by  Sawford  and  Yeung.4,5  These  authors  compared  La¬ 
grangian  autocorrelations  predicted  by  two  dispersion  mod¬ 
els  to  DNS  data.  Both  of  these  models  take  B  to  be  isotropic. 

We  show  here  that  appropriate  values  of  A  and  B  can  be 
deduced  from  the  measured  Lagrangian  velocity  autocovari¬ 
ance,  and  that  the  resulting  model  is  in  good  agreement  with 
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the  DNS  data.  This  agreement  supports  the  nontrivial  con¬ 
clusion  that  the  Lagrangian  velocity  is  well  represented  by  a 
linear  diffusion  process  (except  over  small  time  intervals). 
The  deduced  value  of  B  is  significantly  anisotropic. 

The  performance  of  two  models  used  in  PDF  methods  is 
compared  to  the  DNS  data.  These  are  the  simplified  Lange- 
vin  model  (SLM)  and  the  Lagrangian  isotropization  of  pro- 
duction  (LIPM)  model.14 


II.  HOMOGENEOUS  TURBULENT  SHEAR  FLOW 


In  homogeneous  turbulent  shear  flow,  the  imposed  mean 
velocity  gradient  is 


dJCj 


=  SSiX 


<5/2  > 


(6) 


where  S  is  the  (constant)  imposed  mean  shear  rate.  The  tur¬ 
bulence  is  characterized  by  the  Reynolds  stress  tensor 
( UjUj ),  the  turbulent  kinetic  energy  fc=KM/W/)’  the 
mean  dissipation  rate  s.  All  of  these  quantities  are  uniform  in 
space  and  evolve  in  time. 

An  essential  observation  from  experiments  and  DNS  is 
that,  after  an  initial  transient,  the  turbulence  tends  to  an  ap¬ 
proximately  self-similar  state.  The  normalized  Reynolds- 
stress  tensor 


C 


u 


k 


(7) 


becomes  constant,  as  does  the  ratio  of  turbulence-to-shear 
timescales.  Skis ,  and  hence  also  the  ratio  of  production  V  to 
dissipation  s.  The  turbulent  kinetic  energy  equation  then  dic¬ 
tates  that  k  and  s  increase  exponentially  with  time — as  is 
observed.  Thus  when  normalized  by  k  and  e,  quantities  per¬ 
taining  to  the  energy-containing  scales  of  the  turbulence  are 
self-similar.  Since  the  Reynolds  number  k2/(s  v)  increases 
with  time,  small-scale  quantities  are  not  self-similar  under 
this  scaling. 

The  DNS  of  Sawford  and  Yeung4*5  are  performed  from 
the  nondimensional  time  St—0  until  St  —  20.  The  fluid  par¬ 
ticles  are  introduced  at  <Sr  =  4  when  the  self-similar  state  has 
been  attained.  The  values  Sk/s  =  4.83  and  Vis  — 1.54  are  de¬ 
duced  from  the  values  of  k  and  (u^uj)  from  St  =  4  until  St 
=  20;  and  the  average  value  of  the  normalized  Reynolds 
stress  tensor  over  this  time  interval  is 


C= 


0.96  -0.32  0 

—  0.32  0.43  0 

0  0  0.61 


We  introduce  the  normalized  time 


(8) 


s 


(9) 


and  the  scaled  fluctuating  velocity  following  a  fluid  particle 

do) 

Consistent  with  the  self-similar  state  of  the  turbulence,  we 
assume  that  u (?)  is  a  statistically  stationary  process. 


FIG.  1.  Autocorrelation  functions  p/;(s),  Eq.  (14),  from  the  DNS  data  of 
Sawford  and  Yeung  (Ref.  5). 


The  autocovariance  of  u(0  is 

^ij(s)^(Ui(t)Uj(t  +  s)),  (11) 

which  (in  view  of  the  assumed  stationarity)  is  independent  of 
t ;  and  the  scaled  Reynolds  stress  is 

(ui(t)uj(t))  =  CirRij(0)={-^,  (12) 


which  is  constant.  Note  that  (unlike  Cgj)  Rij(s)  is  not  sym¬ 
metric,  although  it  has  the  property 

4-(*)=£, •,(-*).  (13) 

It  is  conventional  to  define  autocorrelation  functions  by 


Pij(s)= 


Rjjif) 

[C(«)(i)c’0')0)i,/2 


(14) 


(where  bracketed  suffixes  are  excluded  from  the  summation 
convention)  so  that  the  diagonal  components  of  p/;( 0)  are 
unity.  These  autocorrelation  functions  obtained  from  the 
DNS  are  shown  in  Fig.  1.  (Note  that,  by  symmetry,  p23 
=  P  32  =  0-) 

The  analysis  below  shows  that  a  preferable  definition  of 
the  autocorrelations  is 


Rij(s)  =  Cfcl&kj(s),  (15) 

where  C^1  denotes  the  i-k  component  of  the  inverse  of  C. 
Unlike  p/y ,  Rtj  is  a  tensor,  and  at  the  origin  it  is 

Rij(0)=Sij.  (16) 

These  autocorrelation  functions  obtained  from  the  DNS  are 
shown  in  Fig.  2.  [There  is  a  small  inconsistency  in  the  ex¬ 
traction  of  numerical  values  from  the  DNS:  C/y  is  obtained 
as  an  average  from  <5^=4  to  St =20,  whereas  Rgj( 0)  is  ob¬ 
tained  at  St  =  4.  As  a  consequence,  as  may  be  seen  in  Fig.  2, 
the  numerical  values  do  not  satisfy  Eq.  (16)  exactly.] 

Based  on  Rjj(s),  we  define  the  (normalized)  integral 
timescales  by 
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FIG.  2.  Autocorrelation  functions  R^s),  Eq.  (15),  from  the  DNS  data  of 
Sawford  and  Yeung  (Ref.  5). 


r"-J> 


,(s)ds. 


The  values  deduced  from  the  DNS  data  are 


0.44  -  0.06  0 


T=  -0.11  0.22 


III.  STOCHASTIC  MODEL 

The  stochastic  model  considered  is  Eq.  (1)  written  for 

It  is  convenient  to  use  matrix  notation,  and  so  the  equa¬ 
tion  is  written 

-Aud?+BdW,  (19) 

where  (dWdWT)=ldt,  with  I  being  the  identity,  and  T 
denoting  the  transpose. 

The  drift  matrix  A  is  constant  and  it  is  required  that  its 
eigenvalues  have  positive  real  parts.  The  value  of  A  deduced 
from  the  DNS  (below)  has  the  simplest  structure — real  posi¬ 
tive  eigenvalues  and  independent  eigenvectors.  In  this  case  A 
can  be  decomposed  as 

A=VAV-',  (20) 


where  the  columns  of  V  are  the  eigenvectors  of  A,  and  A  is 
the  diagonal  matrix  of  eigenvalues. 

The  diffusion  coefficient  matrix  B  is  also  constant  and, 
without  loss  of  generality,7  we  take  it  to  be  symmetric  (B 
=  Br). 

A.  Autocorrelation  function 

It  is  readily  deduced  from  Eq.  (19)  that  the  autocovari¬ 
ance  matrix  R($)  [Eq.  (11)]  satisfies  the  ordinary  differential 
equation 

dRr 

~j—  =  -A Rr,  for  $2*0.  (21) 


By  post-multiplying  both  sides  of  this  equation  by  C_1,  we 
find  that  R($)  [defined  by  Eq.  (15)]  satisfies  the  same  equa¬ 
tion 

dRT 

—  =  ~AR,  for  $2=0,  *  (22) 

with  the  simple  initial  condition  Rr(0)=I.  The  solution  to 
this  equation  (satisfying  the  initial  condition)  is15 

°°  (—1)” 

Rr($)=exp(-A$)=2)  - 7— A"$",  for  $2^0, 

H= o  «! 


as  may  be  verified  by  differentiating  with  respect  to  $.  It  has 
been  assumed  that  the  eigenvalues  of  A  have  positive  real 
parts,  which  is  a  sufficient  condition  for  exp(-Aj)  to  con- 
verge  to  zero  as  s  tends  to  infinity. 

In  the  case  that  A  has  li nearly  independent  eigenvectors 
the  solution  can  be  written 

Rr($)  =  V exp(  -  A$)  V“ 1 ,  for  $2=0,  (24) 

and  similarly  for  R 

R($)r=Vexp(  — A$)V_1C,  for  $5=0.  (25) 

Thus  each  component  of  the  autocovariance  is  a  linear  com¬ 
bination  of  three  decaying  exponentials — decaying  because 
the  eigenvalues  are  required  to  be  positive. 

For  the  autocorrelation  timescales,  T  [Eq.  (17)]  we  ob- 

tain 

f  oo  r  oo 

Tr=  Jo  Rr(s)ds= Jo  exp(— A$)*/$=A~1.  (26) 

The  conclusion  from  this  development  is  that  the  matrix 
of  autocorrelation  timescales  T  of  the  process  u(f)  generated 
by  the  stochastic  model  Eq.  (19)  is  uniquely  determined  by 
the  drift  matrix  A  as 

T=(A-')r.  (27) 

This  conclusion  depends  on  the  eigenvalues  of  A  having 
positive  real  parts. 

B.  Covariance 

It  follows  from  Eq.  (19)  that  the  covariance  C=(uur) 
evolves  by 

dC 

^-  =  -AC-CAr+BBr.  (28) 

Given  that  B  is  symmetric  and  that  the  process  is  stationary, 
this  leads  to  the  relation 

B2=AC+CAr.  (29) 

C.  Specification  of  stochastic  model  coefficients 

Can  the  model  coefficients  A  and  B  be  chosen  so  that  the 
autocovariance  R($)  from  the  model  matches  that  obtained 
from  DNS  of  homogeneous  turbulent  shear  flow?  Clearly  the 
answer  is  “no,”  since  the  empirical  autocovariances  will  not 
be  of  the  simple  form  implied  by  the  model— i.e.,  sums  of 
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FIG.  3.  Comparison  of  autocorrelation  functions  p,/s),  Eq.  (14),  from  the 
DNS  data  (symbols)  and  from  the  stochastic  model  (lines)  with  coefficients 
determined  from  the  data  [Eq.  (32)  and  Eq.  (33)].  (p22  circles  and  solid  line; 
Pj3  squares  and  dashed  line.) 


three  exponentials.  Nevertheless,  the  preceding  analysis 
shows  that  A  and  B  can  be  chosen  to  match  the  covariance  C 
and  the  timescales  T.  Specifically,  given  T,  A  is  determined 
by 

A=(T-')r  (30) 

[see  Eq.  (27)];  then  B  is  determined  as  the  symmetric  square 
root  of 

B2  =  AC+CAr  (31) 


[see  Eq.  (29)].  Evidently  this  specification  requires  that  T  be 
non  singular.  An  additional  requirement  is  that  T  and  C  be 
such  that  B2  given  by  Eq.  (31)  is  positive  semi-definite. 

For  the  values  of  C  and  T  obtained  from  the  DNS  of 
homogeneous  turbulent  shear  flow,  the  values  of  A  and  B2 
obtained  from  Eq.  (30)  and  Eq.  (31)  are 


A= 


2.45  1.24  0 

0.65  4.90  0 

0  0  4.22 


(32) 


and 


B2  = 


3.90 

-1.18 

0 


— 1.18  0 
3.84  0 

0  5.14 


(33) 


D.  Comparison  of  autocorrelation  functions 

Figure  3  shows  the  comparison  between  the  autocorrela¬ 
tion  functions  Pij(s)  obtained  from  DNS  compared  to  those 
from  the  model  [with  coefficients  given  by  Eq.  (32)  and  Eq. 
(33)].  Inevitably  there  are  qualitative  differences  at  the  ori¬ 
gin,  For  p  j  j ,  for  example,  the  DNS  value  departs  from  unity 
at  the  origin  as  1  ~  pn(s)~s2,  whereas  the  model  departs  as 
l~Pn(s)~M-  This  leads  to  the  model  values  of  Pn(5) 
being  below  the  DNS  values  at  small  times;  and  then,  from 
the  matching  of  the  integral  timescales,  it  is  not  surprising 


that  at  some  later  times  the  model  value  exceeds  the  DNS 
value.  Given  these  inevitable  differences,  the  agreement  be¬ 
tween  the  model  and  the  DNS  is  as  good  as  could  be  ex¬ 
pected.  In  particular  the  model  captures  the  difference  be¬ 
tween  p  i  i  and  the  other  two  diagonal  components  (which  are 
nearly  equal);  and  the  differences  between  p12  and  p21. 


IV.  GENERALIZED  LANGEVIN  MODEL 


In  PDF  methods,  the  stochastic  Lagrangian  model  for 
velocity  that  is  employed  is  the  generalized  Langevin  model 
(GLM).11’12,7  Applied  to  homogeneous  turbulence,  the  model 
for  u(/)  is 

dui=-~^-Ujdt+GijUjdt+(C0e)i,2dWi,  (34) 

where  the  constant  C0  is  generally  ascribed  the  value  2.1. 
The  coefficient  Gi;  can  depend  on  («,«,},  e  and  d{U)ldXj : 
two  particular  specifications  of  Gy  are  considered  below. 

The  transformation  of  Eq.  (34)  to  an  SDE  for  u(f)  re¬ 
sults  in  the  general  stochastic  model,  Eq.  (19),  with  coeffi¬ 
cients 


<5,7  + 


k  d(  Uj) 

S  dXj 


(35) 


and 


«rco%- 


(36) 


Equation  (35)  can  be  rearranged  to  yield  the  value  of 
(kle)Gij  implied  by  the  DNS: 


3.59  0 

—  4.63  0 

0  -3.95 


(37) 


Since  B  is  found  to  be  anisotropic — as  discussed  further  in 
the  next  subsection — no  choice  of  C0  in  Eq.  (36)  yields  the 
correct  diffusion  coefficient.  Nevertheless,  the  magnitude  of 
the  diffusion  is  characterized  by 


C0^  j  trace(B2),  (38) 

the  value  of  which  deduced  from  the  DNS  is  C0  =  4.3.  By 
comparison,  the  standard  model  Eq.  (36)  yields  C0=C0 
=  2.1. 


A.  Anisotropy  of  the  diffusion  coefficient 

The  GLM,  and  also  dispersion  models,  take  the  diffusion 
coefficient  B  to  be  isotropic,  Eq.  (36).  The  reason  generally 
advanced  for  this  specification  is  consistency  with  the  Kol¬ 
mogorov  hypotheses.  For  (dimensional)  time  intervals  s  in 
the  inertial  subrange,  r^s^kfe  (where  tv  is  the  Kolmog¬ 
orov  timescale),  the  Kolmogorov  hypotheses  predict  that  the 
second-order  Lagrangian  structure  function  is  isotropic  and 
linear  in  s ,  i.e., 

([ui(t+s)-ui(t)][uj(t+s)-Uj(t)])=C0ssSij,  (39) 

where  C0  is  a  Kolmogorov  constant.  The  GLM  yields  pre¬ 
cisely  this  results  if  Cq  is  taken  to  be  Cq  . 


1700  Phys.  Fluids,  Vol.  14,  No.  5,  May  2002 


Stephen  B.  Pope 


FIG.  4.  The  scaled  tensors  B2  and  C  shows  as  ellipses  in  the  -x2  plane. 
The  dot-dashed  line  is  jc^yCy2=(C«)“2t  and  the  solid  line  is  the  corre¬ 
sponding  ellipse  for  B2.  Shown  for  reference  are  dashed  lines  at  0°,  45°, 
90°,  and  135°. 


However,  the  value  of  B2  deduced  from  the  DNS  is  de¬ 
cidedly  anisotropic:  the  eigenvalues  of  B2  (which  are  all 
equal  to  C0  =  2.1  in  the  GLM)  are  found  to  be  2.69,  5.06,  and 
5.14.  It  is  possible  that  this  anisotropy  is  a  Reynolds-number 
effect,  which  vanishes  at  sufficiently  high  Reynolds  number. 
This  possibility  could  be  investigated  through  DNS  at  differ¬ 
ent  Reynolds  numbers. 

It  is  also  possible  that  the  anisotropy  in  the  deduced 
value  of  B2  persists  at  high  Reynolds  numbers,  not  because 
the  Kolmogorov  hypothesis  [Eq.  (39)]  is  incorrect,  but  be¬ 
cause  the  stochastic  Lagrangian  model,  Eq.  (19),  is  too 
simple  to  represent  the  multi-timescale  aspects  of  anisotropic 
turbulence. 

Given  the  observation  that  B2  is  anisotropic,  it  is  natural 
to  consider  modifications  to  the  GLM  to  incorporate  such 
anisotropy.  The  natural  way  to  introduce  anisotropy  in  the 
model  is  to  make  the  diffusion  coefficient  dependent  on  the 
normalized  Reynolds  stresses  C.  But  any  such  model  implies 
that  the  principal  axes  of  B2  and  C  are  aligned,  which  is  not 
supported  by  the  data.  Figure  4  shows  the  ellipses  in  the  x{ 
~x2  plane  corresponding  to  the  tensors  B2  and  C.  The  mis¬ 
alignment  of  the  principal  axes  is  evident.  In  fact,  to  within 
1°,  the  minor  axis  of  B2  is  aligned  with  the  major  axis  of  the 
mean  rate-of-stain  tensor  S  (i.e.,  the  45°  line  x2=xl).  Hence 
an  anisotropic  model  for  B2  could  be  constructed  based  on  S 
that  is  consistent  with  the  DNS  data.  More  data — from  dif¬ 
ferent  flows  and  at  different  Reynolds  numbers — are  needed 
before  an  anisotropic  model  for  B2  of  any  generality  can  be 
constructed. 

B.  Simplified  Langevin  model 

In  this  and  the  next  subsection  we  examine  two  specific 
forms  of  the  generalized  Langevin  model,  corresponding  to 
particular  specifications  of  G,y . 

For  the  simplified  Langevin  model  (SLM)  considered 
here,  the  specification  is 


TABLE  I.  Values  of  the  mean  timescale  7=  y  trace(T),  the  normalized 
Reynolds  stress  C;j ,  and  the  turbulence-to-shear  timescale  ratio  Skle  from 
the  DNS  of  Sawford  and  Yeung  (Ref.  5)  and  from  SLM  and  LEPM  for 
different  values  of  C0  and  a2 . 


DNS 

SLM 

SLM 

LIPM 

LIPM 

Q 

- 

2.1 

3.4 

2.1 

4.4 

a2 

- 

- 

- 

3.5 

11.9 

T 

0.30 

0.43 

0.30 

0.63 

0.30 

c„ 

0.96 

1.10 

0.98 

1.02 

1.02 

C22 

0.43 

0.45 

0.51 

0.49 

0.49 

C33 

0.61 

0.51 

0.51 

0.49 

0.49 

^12 

—0.32 

-0.39 

-0.34 

-0.36 

-0.36 

Skle 

4.83 

4.02 

4.47 

4.28 

4.28 

Gij- 


13  \s 
2  +  4  C»)*V 


so  that  the  matrix  A  [Eq.  (35)]  is 


A= 


X  <t  0 

0X0 
0  0  X 


with 


1  V  3 

X"27  +  4C°’ 


(40) 


(41) 


(42) 


and 


a— Skis. 


Evidently  all  three  eigenvalues  of  A  are  equal  to  \,  and  the 
eigenvectors  are  not  independent — two  are  equal  to  [  1  0  0]r. 
Consequently,  the  autocovariance  R (s)  is  not  given  by  Eq. 
(25),  but  instead  the  solution  to  Eq.  (21)  is 


R  (s)  =  e~Xs 


C  ii  crsCi2 
C2{~  crsC22 


0 


C12 

^22 

0 


(43) 


Given  a  specified  value  of  C0  and  the  DNS  value  of 
VI  e,  Eq.  (31)  can  be  solved  to  determine  the  normalized 
Reynolds  stresses  C  given  by  SLM  in  homogeneous  turbu¬ 
lent  shear  flow,  and  then  the  autocovariances  can  be  evalu¬ 
ated  from  Eq.  (43).  We  consider  two  values  of  C0:  the  stan¬ 
dard  value  Cp= 2.1;  and  the  value  Co =3.4  for  which  the 
SLM  timescale  \~l  matches  the  average  DNS  timescale  T 
=  y  trace  (T).  The  values  of  C  obtained  are  shown  in  Table  I. 

The  autocorrelations  pzy(s)  obtained  with  C0  =  3.4  are 
compared  to  the  DNS  data  in  Fig.  5.  As  expected,  the  agree¬ 
ment  is  much  better  with  the  timescales  matched  (C0  =  3.4) 
than  otherwise  (C0  =  2.1,  not  shown).  The  model  correctly 
predicts  the  equality  of  p 22  and  p33  and  their  distinction  from 
Pn  ,  but  the  quantitative  agreement  is  noticeably  poorer  than 
in  Fig.  3. 

The  model  predicts  a  more  substantial  difference  be¬ 
tween  p\2(s)  and  p2l(s)  than  is  evident  in  the  DNS — a  be¬ 
havior  which  is  easily  understood.  The  only  off-diagonal 
term  [a  in  Eq.  (41)]  enters  the  SDE  for  velocity  as 
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FIG.  5.  Comparison  of  autocorrelation  functions  pl7(s),  Eq.  (14),  from  DNS 
(symbols)  and  from  SLM  with  C0-  3.4.  (For  SLM,  p22=P33  ■) 


FIG.  6.  Comparison  of  autocorrelation  functions  p7(s),  Eq.  (14),  from  DNS 
(symbols)  and  from  LIPM  with  C0=4.4  and  a2=  1 1.9  (lines).  Dashed  lines, 
pl2  and  P33 ;  solid  lines,  p2i  and  p22  • 


dui  =  — — u2dt ....  (44) 

Thus  large  positive  (or  negative)  values  of  u2  tend  to  lead  to 
large  negative  (or  positive)  values  of  ux  after  a  time  lag. 
Thus  the  peak  correlation  \{u2(i)&\(t +•*))! — °r  equiva¬ 
lently  the  minimum  of  p2j(s) — occurs  for  a  positive  value  of 
s.  The  model  overestimates  this  effect,  because  it  takes  no 
account  of  rapid  pressure  fluctuations  which  tend  to  counter¬ 
act  the  effects  of  mean  shear. 


C.  Lagrangian  isotropization  of  production  model 
(LIPM) 

The  LIPM14,7  corresponds  closely  to  the  Launder,  Reece, 
and  Rodi16  Reynolds-stress  model.  Using  standard  values  for 
the  model  constants  /?i~3  and  _6?  the  LIPM  equation  for 
G  is 


-G=a1I+a2(b-3b2) 

e 


Sk 


3 

~5bl2 
1  3 

~5~5*22 


4  3 

5+5b 11 


'12 


(45) 


0  0  i 

where  b  is  the  anisotropy  tensor 

b=ic-^I,  (46) 

2  3 

the  standard  value  of  the  constant  £*2  is  £*2  =  3.5,  and  the 
coefficient  a{  is  given  by 


1 


3  V 


+ jC0  +  — -  +  3a2trace(bJ). 


(47) 


With  the  standard  value  Co  — 2.1,  the  model  yields  rea¬ 
sonable  values  of  the  normalized  Reynolds  stresses,  but  the 
average  time  scale  \  trace (T)  is  more  than  twice  the 


DNS  value;  see  Table  I.  As  a  consequence,  the  model  (with 
C0  =  2.1)  produces  autocorrelations  Pij(s)  (not  shown)  in 
very  poor  agreement  with  the  DNS  data. 

To  provide  a  more  meaningful  comparison,  the  constants 
C0  and  a2  are  adjusted  to  match  the  average  timescale,  while 
leaving  the  normalized  Reynolds  stresses  the  same.  The  au¬ 
tocorrelations  given  by  LIPM  with  these  values  (C0  =  4.4, 
£*2=11.9)  are  compared  to  the  DNS  data  in  Fig.  6.  The 
agreement  is  quite  poor.  Except  at  small  times,  p22{s)  is 
incorrectly  predicted  to  be  larger  than  p33(s);  and  evidently 
the  effect  of  the  rapid  pressure  is  overpredicted  as  there  is 
little  difference  between  pi2(s)  and  p2i(s). 

This  last  point  can  be  seen  directly  in  the  matrix  A, 
which  for  LIPM  is 


A= 


2.86  1.99  0 

2.21  6.11  0 

0  0  4.48 


(48) 


The  direct  effect  of  shear  appears  in  the  1-2  component,  and 
in  SLM  the  2-1  component  is  zero.  In  A  deduced  from  the 
DNS  data  [Eq.  (32)],  A21  as  about  half  of  A12;  but  for  LIPM 
A21  exceeds  A12. 


V.  CONCLUSIONS 

As  previously  observed  by  Sawford  and  Yeung4’5  in  the 
context  of  turbulent  dispersion,  Lagrangian  data  from  DNS 
of  homogeneous  turbulence  is  valuable  in  the  development 
and  testing  of  stochastic  Lagrangian  models.  After  an  initial 
transient,  homogeneous  turbulent  shear  becomes  (approxi¬ 
mately)  self-similar,  so  that  the  appropriately  scaled  La¬ 
grangian  velocity  fluctuation  u (?)  becomes  a  statistically  sta¬ 
tionary  random  process. 

The  stochastic  Lagrangian  model  considered  for  u (t)  is 
the  diffusion  process  Eq.  (19)  in  which  the  drift  coefficient 
depends  linearly  on  u (?)  through  the  drift  matrix  A,  and  the 
(anisotropic)  diffusion  coefficient  B  is  constant.  An  analysis 
of  this  model  shows  that  there  is  a  unique  specification  of  A 
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and  B  [Eq.  (30)  and  Eq.  (31)]  such  that  the  covariance  matrix 
C  and  the  timescale  matrix  T  match  those  obtained  from 
DNS.  The  autocorrelation  functions  predicted  by  the  model 
are  in  good  agreement  with  the  DNS  data  (except  at  small 
times).  The  model  for  u(?)  is  a  continuous,  Gaussian,  Mar¬ 
kov  process;  and  it  is  a  significant  conclusion  that  such  a 
simple  process  provides  a  good  model  for  the  Lagrangian 
velocity  in  homogeneous  turbulent  shear  flow.  (It  is  known 
that  the  one-point  one-time  joint  PDF  of  velocity  is  jointly 
normal1  in  this  flow.) 

Contrary  to  conventional  modelling  assumptions,  if  is 
found  that  the  diffusion  coefficient  B  is  significantly  aniso¬ 
tropic.  Whether  or  not  this  is  a  low  Reynolds-number  effect 
is  an  important  question  which  can  be  addressed  in  future 
DNS  studies. 

The  magnitude  of  the  diffusion  coefficient  can  be  char¬ 
acterized  by  C0— j  trace(B^)  and  the  value  deduced  from 
the  DNS  data  is  C0  =  4.3.  This  is  substantially  larger  than  the 
corresponding  value  C0  =  2.1  normally  used  in  PDF  models. 

There  is  evidence  that  the  appropriate  value  of  C0  de¬ 
pends  on  Reynolds  number.13’10  In  the  DNS,  the  Taylor-scale 
Reynolds  number  based  on  -direction  statistics  increases 
from  RK^40  to  R^llO  during  the  course  of  the  simulation. 
Sawford  and  Yeung5  provide  an  empirical  expression  for  Cq 
as  a  function  of  /?x,  which  increase  from  C0  =  3.7  at  Rk 
“40  to  C0  =  5.4  at  /?x=110.  The  value  C0=4.3  deduced 
from  the  DNS  lies  within  this  range,  but  Reynolds-number 
effects  are  not  addressed  here. 

The  autocorrelation  functions  predicted  by  two  general¬ 
ized  Langevin  models  are  compared  to  the  DNS  data  in  Figs. 
5  and  6.  In  the  simplified  Langevin  model  (SLM),  no  ac¬ 
count  is  taken  of  the  rapid  pressure  fluctuations,  and  as  a 
consequence  the  difference  between  pi2(s)  and  p2i(s)  is 
overpredicted.  The  Lagrangian  IP  model  (LIPM) — which  in¬ 
cludes  a  model  for  the  rapid  pressure — yields  autocorrela¬ 
tions  in  poor  agreement  with  the  DNS  data. 

In  the  specification  of  both  the  drift  and  diffusion  coef¬ 
ficients,  there  is  clearly  scope  for  considerable  improvement 
in  generalized  Langevin  models. 
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A  methodology  termed  the  “velocity  filtered  density  function”  (VFDF)  is  developed  and 
implemented  for  large  eddy  simulation  (LES)  of  turbulent  flows.  In  this  methodology,  the  effects  of 
the  unresolved  subgrid  scales  (SGS)  are  taken  into  account  by  considering  the  joint  probability 
density  function  of  all  of  the  components  of  the  velocity  vector.  An  exact  transport  equation  is 
derived  for  the  VFDF  in  which  the  effects  of  the  SGS  convection  appear  in  closed  form.  The 
unclosed  terms  in  this  transport  equation  are  modeled.  A  system  of  stochastic  differential  equations 
(SDEs)  which  yields  statistically  equivalent  results  to  the  modeled  VFDF  transport  equation  is 
constructed.  These  SDEs  are  solved  numerically  by  a  Lagrangian  Monte  Carlo  procedure  in  which 
the  Ito-Gikhman  character  of  the  SDEs  is  preserved.  The  consistency  of  the  proposed  SDEs  and  the 
convergence  of  the  Monte  Carlo  solution  are  assessed  by  comparison  with  results  obtained  by  an 
Eulerian  LES  procedure  in  which  the  corresponding  transport  equations  for  the  first  two  SGS 
moments  are  solved.  The  VFDF  results  are  compared  with  those  obtained  via  several  existing  SGS 
closures.  These  results  are  also  analyzed  via  a  priori  and  a  posteriori  comparisons  with  results 
obtained  by  direct  numerical  simulation  of  an  incompressible,  three-dimensional,  temporally 
developing  mixing  layer.  ©  2002  American  Institute  of  Physics.  [DOI:  10.1063/1.1436496] 


I.  INTRODUCTION 

The  probability  density  function  (PDF)  approach  has 
proven  useful  for  large  eddy  simulation  (LES)  of  turbulent 
reacting  flows.1-15  The  formal  means  of  conducting  such 
LES  is  by  consideration  of  the  “filtered  density  function” 
(FDF)  which  is  essentially  the  filtered  fine-grained  PDF  of 
the  transport  quantities.  In  all  previous  contributions,  the 
FDF  of  the  “scalar”  quantities  is  considered:  Gao  and 
O’Brien,3  Colucci  et  al.}  Reveillon  and  Vervisch,7  and  Zhou 
and  Pereira12  developed  a  transport  equation  for  the  FDF  in 
constant  density  turbulent  reacting  flows.  Jaberi  et  al}  ex¬ 
tended  the  methodology  for  LES  of  variable  density  flows  by 
consideration  of  the  “filtered  mass  density  function” 
(FMDF),  which  is  essentially  the  mass  weighted  FDR  The 
fundamental  property  of  the  PDF  methods  is  exhibited  by  the 
closed  form  nature  of  the  chemical  source  term  appearing  in 
the  transport  equation  governing  the  FDF  (FMDF).  This 
property  is  very  important  as  evidenced  in  several  applica¬ 
tions  of  FDF  for  LES  of  a  variety  of  turbulent  reacting 
flows.6-10,12  However,  since  the  FDF  of  only  the  scalar  quan¬ 
tities  are  considered,  all  of  the  “hydrodynamic”  effects  are 
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modeled.  In  all  previous  LES/FDF  simulations,  these  effects 
have  been  modeled  via  “non-FDF”  methods. 

The  objective  of  the  present  work  is  to  extend  the  FDF 
methodology  to  also  include  the  subgrid  scale  (SGS)  velocity 
vector.  This  is  facilitated  by  consideration  of  the  joint  “ve¬ 
locity  filtered  density  function”  (VFDF).  With  the  definition 
of  the  VFDF,  the  mathematical  framework  for  its  implemen¬ 
tation  in  LES  is  established.  A  transport  equation  is  devel¬ 
oped  for  the  VFDF  in  which  the  effects  of  SGS  convection 
are  shown  to  appear  in  closed  form.  The  unclosed  terms  in 
this  equation  are  modeled  in  a  fashion  similar  to  that  in  the 
Reynolds-averaged  simulation  (RAS)  procedures.  A  La¬ 
grangian  Monte  Carlo  procedure  is  developed  and  imple¬ 
mented  for  numerical  simulation  of  the  modeled  VFDF 
transport  equation.  The  consistency  of  this  procedure  is  as¬ 
sessed  by  comparing  the  first  two  moments  of  the  VFDF 
with  those  obtained  by  the  Eulerian  finite  difference  solu¬ 
tions  of  the  same  moments  transport  equations.  The  results  of 
the  VFDF  simulations  are  compared  with  those  predicted  by 
the  Smagorinsky16  closure,  and  the  “dynamic”  Smagorinsky 
model.17-19  The  VFDF  results  are  also  assessed  via  compari¬ 
sons  with  direct  numerical  simulation  (DNS)  data  of  a  three- 
dimensional  (3D)  temporally  developing  mixing  layer  in  a 
context  similar  to  that  of  Vreman  et  al 20 

This  work  deals  with  LES  of  the  velocity  field  in  a  con¬ 
stant  density,  nonreacting  flow.  Consideration  of  the  joint 
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Note  that  these  models  (i.e.,  the  first  two  terms  on  the  right- 
hand  sides  of  Eqs.  (18)  and  (19)  are  the  same,  but  that  they 
model  slightly  different  quantities.  With  this  closure,  the  two 
terms  in  GI;  and  e  jointly  represent  the  SGS  pressure-strain 
and  SGS  dissipation.  These  are  modeled  as11,26 

G,7=-£o(^+^C°)<5,7,  E  =  Csky2/AL,  o)  =  e/k, 

(20) 


where  cj  is  the  SGS  mixing  frequency,  A,  is  the  filter  width, 
k  =  jrL(Uj ,«,)  is  the  SGS  kinetic  energy,  and  e=jeii  is  the 
SGS  dissipation  rate. 

With  the  GLM,  the  two  forms  of  the  VFDF  transport 
equation  are 


DPl 
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d 

dxk 
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i  d*P 

[C, ,(.,-(«, U/>J+rCoe—i:.  (21) 


dVj 

for  VFDF1,  and 


2  "  dvtdVi 

3{p)l  SPl 


djVjklL  dP l 
dxk  dVi 


lGij(Vj-(Uj)L)PL 1 
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for  VFDF2.  Hereinafter,  Eqs.  (21)  and  (22)  are  referred  to  as 
“VFDF1”  and  “VFDF2,”  respectively.  The  difference  be¬ 
tween  these  two  equations  is  in  the  different  treatment  of  the 
closed  viscous  terms. 


D.  Transport  equations  for  moments 

The  zeroth,  first,  and  second  moment  equations  corre¬ 
sponding  to  these  two  formulations  are 

for  VFDF1: 
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for  VFDF2: 
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It  may  be  seen  that  the  zeroth  and  first  moment  equations  are 
identical  (and  exact);  whereas  the  second  central  moment 
equations  differ  by  the  additional  viscous  term  in  VFDF1 
[Eq.  (24)].  A  comparison  of  these  modeled  equations  with 
Eq.  (5)  shows  that  the  GLM  model  implies 

-  FI  ij  -  ( 8  ir  f e  <?,7)  =  -  Cx  <•>[  tl(  u  i ,  Uj)  - 1 k  8ij ] , 

C!=1  +  |C0.  (27) 

This  is  the  same  as  the  Rotta33  model  as  shown  by  Pope.34 
There  are  two  model  constants  in  the  VFDF  equation.  In 
RAS,  typically34,35  Ce«*l,  and  C0~2.1  (Ci  =  4.15).  As 
shown  in  Refs.  27,  34  boundedness  of  the  GLM  coefficients 
C0> 0  guarantees  that  the  SGS  stress  is  realizable. 


IV.  EQUIVALENT  STOCHASTIC  SYSTEMS 

The  solution  of  the  VFDF  transport  equation  provides  all 
the  statistical  information  pertaining  to  the  velocity  vector. 
The  most  convenient  means  of  solving  this  equation  is  via 
the  Lagrangian  Monte  Carlo  scheme.  The  basis  of  this 
scheme  relies  upon  the  principle  of  equivalent  systems.26,32 
Two  systems  with  different  instantaneous  behaviors  may 
have  identical  statistics  and  satisfy  the  same  PDF  transport 
equation.  In  this  context,  the  general  diffusion  process  is 
considered  via  the  following  system  of  stochastic  differential 
equations  (SDEs):26’31’36'37 

dX{t)  =  Di{Z{t)M(.t)\t)dt  +  B(X(t),U{t)\t)dWxi(t), 

dUi(t)  =  Mi(X(t)M(t);t)dt+E(X(t)Mt);t)dWvi(t) 

+  Fij(X(t)M(t);t)dWx(t),  (28) 

where  and  ZY,  are  probabilistic  representations  of  x  and  w, 
respectively.  The  coefficients  Dt  and  Mt  are  the  “drift”  in 
the  phase  space  of  position  and  velocity,  respectively.  The 
terms  B  and  E  are  the  “diffusion”  coefficients  for  physical 
and  velocity  spaces,  respectively;  and  W *  and  W”  denote 
independent  Wiener-Levy  processes.38  The  tensor  F \j  repre¬ 
sents  the  dependency  between  the  velocity  and  physical 
spaces.  This  term  is  needed  to  satisfy  the  Ito  condition  for 
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B^O.  A  comparison  of  the  Fokker-Planck  equation  of  Eq. 
(28)  with  the  modeled  VFDF1  transport  equation,  Eq.  (21) 
yields 


d{ p)L  d2(ui)L 

M‘—-^+2v^rk+G‘^rW'  °<-U. 


B  —  yj2v ,  E—  -\JCqe  ,  F ij—  \l2  v* 


d{Ui) 


(29) 


i/L 


Bx  ; 


Therefore,  the  proper  SDEs  which  represent  VFDF1  in  the 
Lagrangian  sense  are 
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This  stochastic  system  is  the  same  as  that  developed  by 
Dreeben  and  Pope29-31  for  RAS. 

For  VFDF2,  due  to  the  absence  of  diffusion  in  physical 
space  we  must  have  B-  0.  Therefore,  the  corresponding 
SDEs  are 


dx^uMdt, 


dUi{t)  = 
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d{°jkh. 

dxk 
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dt 
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This  system  is  the  same  as  that  suggested  by  Pope26  and 
Haworth  and  Pope27  for  RAS. 

The  primary  difference  between  the  two  formulations 
VFDF1  and  VFDF2  is  due  to  molecular  effects  in  the  spatial 
diffusion  of  the  VFDF.  This  is  explicitly  included  in  the 
VFDF1  formulation  and  is  also  present  in  the  corresponding 
second  moment  equation.  This  difference  is  expected  to  be 
important  in  flows  where  viscous  effects  are  important;  e.g., 
flow  near  solid  boundaries.29-31  Both  of  these  formulation 
are  considered  in  our  numerical  simulations  as  discussed  be¬ 
low. 


V.  NUMERICAL  SOLUTION  PROCEDURE 

Numerical  solution  of  the  modeled  VFDF  transport 
equation  is  obtained  by  a  Lagrangian  Monte  Carlo  proce¬ 
dure.  The  basis  of  this  procedure  is  the  same  as  that  in 
RAS39-41  and  in  previous  LES/FDF.6,8  But  there  are  some 
subtle  differences  which  are  explained  here.  In  the  Lagrang¬ 
ian  description,  the  VFDF  is  represented  by  an  ensemble  of 
N  statistically  identical  Monte  Carlo  particles.  Each  of  these 
particles  carries  information  pertaining  to  its  velocity 
U*n*(0  and  position  X^(/),  n  —  1,2, ...N.  This  information  is 
updated  via  temporal  integration  of  Eq.  (28).  The  simplest 
means  of  performing  this  integration  is  via  the  Euler- 
Maruyamma  approximation42 


X?(/Jk+1)  =  X?(^)  +  D”(/,)  ±t  +  B»(tk)(At)mg(tk), 

f/?(^+i)  =  ^(4)+A/?(tt)Ar+£"(/t)(A01/2^(/t) 

+  Flj(tk)(At)mCl(tk),  (32) 

where  D?(tk)  =  Di(X<n\tk),V<n\tky,t),  B<">(r*) 
-B(X<-n)(tk),lJM(tk);t),...  and  i;?(tk),  £](tk)  are  indepen- 
dent  standardized  Gaussian  random  variables.  This  formula¬ 
tion  preserves  the  Markovian  character  of  the  diffusion 
processes43,44  and  facilitates  affordable  computations. 
Higher-order  numerical  schemes  for  solving  Eq.  (28)  are 
available,42  but  one  must  be  cautious  in  using  them  for  LES.6 
Since  the  diffusion  term  in  Eq.  (28)  strongly  depends  on  the 
stochastic  processes,  the  numerical  scheme  must  be  consis¬ 
tent  with  Ito-Gikhman45,46  calculus.  Equation  (32)  exhibits 
this  property. 

The  statistics  are  evaluated  by  consideration  of  the  en¬ 
semble  of  particles  in  a  “finite  volume”  centered  at  a  spatial 
location.  This  ensemble  provides  “one-time”  statistics.  This 
finite  volume  is  characterized  by  a  cubic  box  of  length  A£ . 
This  is  necessary  as,  with  probability  one,  no  particle  will 
coincide  with  the  point  as  considered.32  Here,  a  cubic  box  of 
size  A£  is  used  to  construct  the  ensemble  mean,  variances 
and  covariances  of  the  velocity  vector.  These  values  are  used 
in  the  finite  difference  LES  solver  of  Eq.  (4)  as  described 
below. 

The  SGS  dissipation  rate  and  the  SGS  mixing  frequency 
as  required  in  the  solution  of  the  VFDF  are  evaluated  on  the 
finite  difference  grid  points  and  interpolated  to  the  particle’s 
location.  Ideally,  for  reliable  Eulerian  statistics  and  minimum 
numerical  dispersion,  it  is  desired  to  have  the  size  of  the 
sample  domain  infinitesimally  small  (i.e.,  A£—>0)  and  the 
number  of  particles  within  this  domain  infinitely  large.  That 
is 


PL(v;x,t)  + - VN  (v;x,t)=2-  X  £(i>-w(n)), 

N  E  ne AE 


Nf— 

Ae~>0 


(33) 


where  VNe  is  the  Eulerian  PDF  constructed  from  the  particle 
ensemble,  neA£  denotes  the  particles  contained  in  an  en¬ 
semble  box  of  length  A£  centered  at  x;  and  NE  is  the  total 
number  of  particles  within  the  box.  With  a  finite  number  of 
particles,  obviously  a  larger  A£  is  needed.  This  compromise 
between  the  statistical  accuracy  and  dispersive  accuracy  im¬ 
plies  that  the  optimum  magnitude  of  A£  cannot,  in  general, 
be  specified  a  priori. 11,26  This  does  not  diminish  the  capabil¬ 
ity  of  the  procedure,  but  exemplifies  the  importance  of  the 
parameters  governing  the  statistics. 

To  provide  an  estimate  of  the  proper  A£  size,  a  “point 
estimator”  procedure  is  considered.  With  this  procedure,  the 
mean  values  (the  first  moments  of  the  VFDF)  are  evaluated 
by  ensemble  averaging,  and  spatial  variations  of  these  mean 
values  within  the  box  are  ignored.  With  the  discrete  repre¬ 
sentation  [Eq.  (32)],  the  first  two  moments  in  this  procedure 
are  evaluated  via 
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TABLE  I.  Recapitulation  of  the  VFDF  solution  procedures. 


Finite 

Particle 

Particle  statistics 

Finite  difference 

difference 

solver 

used  by  the  finite 

variables  used  by 

Redundant 

variables 

variables 

difference  solver 

particle  solver 

quantities 

VFDF  1 

Ml 

x, 

rL(uitUj) 

. ,  <Kp>i 

Ml 

{ p)l 

Vi 

tot  Auk 

dxk  ’  dxkdxk 

VFDF  2 

Ml 

X, 

rL(Ui,Uj ) 

(u,)l 

(p)l 

U, 

Auk 

dxkdxk 

LES-FD 

<“. )lAp)l 

X, 

TL(uiyUj,Uk) 

/  s  A>)l 

Ml 

Tl(ui  >“/) 

Vi 

Auk 

Tl(Ui,Uj) 

dxkdxk 

(“<>L 


Nr--*00 

0 


£ * 


TL(Ui,Uj) 


Nr—*00 

ae-*o 


(^  -  <  c/,-)  ( r/j"^  -  <  t/y)£) . 

(34) 


The  point  estimator  is  obviously  subject  to  both  statistical 
errors  and  dispersive  errors  for  A^O. 

To  determine  the  pressure  field,  the  “mean-field  solver” 
is  based  on  the  “compact  parameter”  finite  difference 
scheme  of  Carpenter.47  This  is  a  variant  of  the  McCormack48 
scheme  in  which  fourth-order  compact  differences  are  used 
to  approximate  the  spatial  derivatives,  and  a  second-order 
symmetric  predictor-corrector  sequence  is  employed  for  time 
discretization.  The  numerical  algorithm  is  a  hyperbolic 
solver  which  considers  a  fully  compressible  flow.  Here,  the 
simulations  are  conducted  with  a  low  Mach  number  (M 
^0.3)  to  minimize  compressibility  effects.  All  the  finite  dif¬ 
ference  operations  are  conducted  on  fixed  and  equally  sized 
grid  points.  The  transfer  of  information  from  these  points  to 
the  location  of  the  Lagrangian  particles  is  conducted  via  in¬ 
terpolation.  A  second-order  (bilinear)  interpolation  scheme  is 
used  for  this  purpose.  The  results  of  previous  work  indicate 
no  significant  improvements  with  the  use  of  higher  order 
interpolation  schemes.6 

The  mean-field  solver  also  determines  the  filtered  veloc¬ 
ity  field.  That  is,  there  is  a  “redundancy”  in  the  determina¬ 
tion  of  the  first  filtered  moments  as  both  the  finite  difference 
and  the  Monte  Carlo  procedures  provides  the  solution  of  this 
field.  This  redundancy  is  actually  very  useful  in  monitoring 
the  accuracy  of  the  simulated  results.  Detailed  discussions 
pertaining  to  this  issue  are  provided  in  Refs.  8,  39-41. 

To  establish  the  consistency  of  the  VFDF  solver,  another 
LES  is  also  conducted  in  which  the  modeled  transport  equa¬ 
tions  for  the  filtered  velocity  and  the  generalized  SGS 
stresses  are  solved  strictly  via  the  finite  difference  scheme. 
These  simulations  are  referred  to  as  LES-FD  and  are  only 


applied  for  the  case  corresponding  to  VFDF2.  That  is,  Eqs. 
(25)  and  (26)  are  considered.  Since  the  SGS  transport  terms 
rL{Ui,Uj,uk)  are  unclosed  in  Eq.  (26),  the  values  corre¬ 
sponding  to  these  terms  are  taken  from  the  Monte  Carlo 
solver  and  substituted  in  the  SGS  stress  transport  equations. 
The  attributes  of  all  of  the  scheme  are  summarized  in  Table  I, 
with  further  discussions  in  Refs.  6,  39-41. 

VI.  RESULTS 
A.  Flows  simulated 

Simulations  are  conducted  of  a  two-dimensional  (2D) 
planar  jet,  and  a  3D  temporally  developing  mixing  layer.  The 
jet  flow  simulations  are  conducted  primarily  for  establishing 
the  consistency  of  the  Lagrangian  Monte  Carlo  solver.  For 
this  purpose,  2D  simulations  are  sufficient.  To  analyze  the 
overall  performance  of  the  VFDF  and  to  demonstrate  its  full 
capabilities  and  drawbacks,  3D  simulations  are  required. 

In  the  planar  jet,  a  fluid  issues  from  a  jet  of  width  D  into 
a  co-flowing  stream  with  a  lower  velocity.  The  size  of  the 
domain  in  the  streamwise  (x)  and  cross-stream  (y)  directions 
are  0=^x=S  14D  and  -3.5D^y^3.5D.  The  ratio  of  the  co¬ 
flowing  stream  velocity  to  that  of  the  jet  at  the  inlet  is  kept 
fixed  at  0.5.  A  double-hyperbolic  tangent  profile  is  utilized  to 
assign  the  velocity  distribution  at  the  inlet  plane.  The  forma¬ 
tion  of  the  large  scale  coherent  structures  are  expedited  by 
imposing  low  amplitude  perturbations  at  the  inlet.  In  the  fi¬ 
nite  difference  simulations,  the  characteristic  boundary  con¬ 
dition  procedure  of  Ref.  49  is  used  at  the  inlet,  free-shear 
boundary  conditions  are  used  at  the  free-streams  and  the 
pressure  boundary  condition  of  Ref.  50  is  used  at  the  outflow. 

The  temporal  mixing  layer  consists  of  two  parallel 
streams  traveling  in  opposite  directions  with  the  same 
speed.51-53  A  hyperbolic  tangent  profile  is  utilized  to  assign 
the  velocity  distribution  at  the  initial  time.  The  simulations 
are  conducted  for  a  cubic  box,  O^x^L,  —  L/2^y=^L/2,  0 
where  x,  y,  and  z  denote  the  streamwise,  the  cross¬ 
stream  and  the  spanwise  directions,  respectively;  and  the 
length,  L  is  specified  such  that  L  —  2Np\u,  where  NP  is  the 
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desired  number  of  successive  vortex  pairings  and  XM  is  the 
wavelength  of  the  most  unstable  mode  corresponding  to  the 
mean  streamwise  velocity  profile  imposed  at  the  initial 
time.  The  flowfield  is  parameterized  in  a  procedure  some¬ 
what  similar  to  that  by  Vreman  et  al 20  The  formation  of  the 
large-scale  structures  are  expedited  through  eigenfunction 
based  initial  perturbations.54,55  This  includes  two- 
dimensional20,52,56  and  three-dimensional52,57  perturbations 
with  a  random  phase  shift  between  the  3D  modes.  This  re¬ 
sults  in  the  formation  of  two  successive  vortex  pairings  and 
strong  three-dimensionality. 

The  flow  variables  are  normalized  with  respect  to  se¬ 
lected  reference  quantities.  In  the  jet  flow,  the  jet  exit  veloc¬ 
ity,  and  the  jet  width  are  the  reference  scales.  In  the  temporal 
mixing  layer,  the  reference  length  is  the  half  initial  vorticity 
thickness,  Lr=Sv(t~ 0)12  (Sv  —  AUlld^^J  dyl^,  where 
{u\)i  is  the  Reynolds  averaged  value  of  the  filtered  stream- 
wise  velocity  and  A  U  is  the  velocity  difference  across  the 
layer).  The  reference  velocity  is  Ur—  AUH. 

B.  Numerical  specifications 

All  finite  difference  simulations  are  conducted  on 
equally  spaced  grid  points  with  grid  spacings  Ax  — Ay 
~Az  ( for  3D)  —  A .  The  resolution  for  LES  of  the  planar  jet 
consists  of  201X101  grid  points.  This  allows  simulations 
with  a  Reynolds  number  Re  =  U rD!v  —  14,000.  The  simula¬ 
tions  of  the  temporal  mixing  layer  are  conducted  on  1933  and 
333  points  for  DNS  and  LES,  respectively.  This  allows  simu¬ 
lations  with  Re  —  UrLr/v  —  50. 

To  filter  the  DNS  data,  a  top-hat  function21  of  the  form 
below  is  used 


G(x'-x)  =  Il  G(x;-x,), 

i-  1 


G(x'i-xi)  = 


A  l*(  *il 


(35) 


A/ 


in  which  ND  denotes  the  number  of  dimensions,  and  AL 
—  2 A. 58  No  attempt  is  made  to  investigate  the  sensitivity  of 
the  results  to  the  filter  function24  or  the  size  of  the  filter.59 

For  VFDF  simulations  of  the  temporal  mixing  layer,  the 
Monte  Carlo  particles  are  initially  distributed  throughout  the 
computational  region.  For  the  jet  flow,  the  particles  are  sup¬ 
plied  in  the  inlet  region  -  1.75D^y^  1.75D.  As  the  par¬ 
ticles  convect  downstream,  this  zone  distorts  as  it  conforms 
to  the  flow  as  determined  by  the  hydrodynamic  field.  The 
simulation  results  are  monitored  to  ensure  the  particles  fully 
encompass  and  extend  well  beyond  regions  of  nonzero  vor¬ 
ticity  with  an  approximately  uniform  particle  number  den¬ 
sity.  All  simulations  are  performed  with  a  uniform 
“weight”26  of  the  Monte  Carlo  particles.  In  the  temporal 
mixing  layer,  due  to  flow  periodicity  in  the  streamwise  and 
spanwise  directions,  if  the  particle  leaves  the  domain  at  one 
of  these  boundaries  new  particles  are  introduced  at  the  other 
boundary  with  the  same  compositional  values.  In  the  cross¬ 


stream  directions,  the  free-slip  boundary  condition  is  satis¬ 
fied  by  the  mirror-reflection  of  the  particles  leaving  through 
these  boundaries.  In  the  planar  jet,  new  particles  are  intro¬ 
duced  through  the  inlet  boundary  at  a  rate  proportional  to  the 
local  flow  velocity  and  with  a  velocity  makeup  dependent  on 
the  cross-stream  direction  only.  When  the  particles  leave  the 
computational  domain  at  the  outflow,  they  are  no  longer 
tracked.  The  density  of  the  Monte  Carlo  particles  is  deter¬ 
mined  by  the  average  number  of  particles  NE  within  the 
ensemble  domain  of  size  A£XA£(XA£).  The  effects  of 
both  of  these  parameters  are  assessed  to  ensure  the  consis¬ 
tency  and  the  statistical  accuracy  of  the  VFDF  simulations. 

All  results  are  analyzed  both  “instantaneously”  and 
“statistically.”  In  the  former,  the  instantaneous  contours 
(snap-shots)  and  scatter  plots  of  the  variables  of  interest  are 
analyzed.  In  the  latter,  the  “Reynolds-averaged”  statistics 
constructed  from  the  instantaneous  data  are  considered.  In 
the  spatially  developing  flows  this  averaging  procedure  is 
conducted  via  sampling  in  time.  In  the  temporal  mixing 
layer,  the  statistics  are  constructed  by  spatial  averaging  over 
the  x-z  plane  of  statistical  homogeneity.  All  Reynolds  aver¬ 
aged  results  are  denoted  by  an  overbar. 

C.  Consistency  and  convergence  assessments 

The  objective  of  this  section  is  to  demonstrate  the  con¬ 
sistency  of  the  VFDF  formulation  and  the  convergence  of  its 
Monte  Carlo  simulation  procedure.  For  this  purpose,  the  re¬ 
sults  via  VFDF  and  LES-FD  are  compared  against  each 
other.  Since  the  accuracy  of  the  finite  difference  procedure  is 
well-established  (at  least  for  the  first-order  filtered  quanti¬ 
ties),  such  a  comparative  assessment  provides  a  good  means 
of  assessing  the  performance  of  the  Monte  Carlo  solution  of 
the  VFDF.  To  do  so,  the  statistical  results  obtained  from  the 
Monte  Carlo  simulations  of  Eq.  (31)  are  compared  with  the 
finite  difference  solution  of  Eqs.  (25)  and  (26).  Also,  no  at¬ 
tempt  is  made  to  determine  the  appropriate  values  of  the 
model  constants;  the  values  suggested  in  the  literature  are 
adopted34  C0=2.1(C1  =  4.15)  and  Ce=  1. 

In  Fig.  1,  the  instantaneous  contour  plots  of  the  vorticity 
are  shown  as  determined  by  (a)  VFDF2  and  (b)  LES-FD. 
This  figure  provides  a  simple  visual  demonstration  of  the 
consistency  of  the  VFDF2.  Scatter  plots  of  (u)L  vs  (v)L  are 
presented  in  Fig.  2.  The  correlation  and  regression  coeffi¬ 
cients  (denoted,  respectively,  by  p  and  r  on  these  figures)  are 
insensitive  to  AE .  Figures  3  and  4  show  the  Reynolds  aver¬ 
aged  values  of  the  streamwise  velocity  and  several  compo¬ 
nents  of  the  SGS  stress  tensor  for  several  values  of  A£ ,  with 
Ne~ 40  kept  fixed.  It  is  observed  that  the  first  filtered  mo¬ 
ments  as  obtained  by  VFDF  agree  very  well  with  those  via 
LES-FD  even  for  large  A£  values.  However,  smaller  A£  val¬ 
ues  are  required  for  convergence  of  the  VFDF  predicted  SGS 
stresses  to  those  by  LES-FD.  The  relative  difference  between 
the  L2  norms  of  all  of  the  components  of  the  SGS  tensor  as 
a  function  of  (A£/A)2  is  presented  in  Fig.  5.  Extrapolation 
to  A£=0  shows  that  the  “error”  goes  to  zero  as  A£— >0. 

The  influence  of  NE  on  the  first  two  moments  is  shown 
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(a) 


FIG.  1.  Plot  of  the  vorticity  field  con¬ 
tours,  (a)  VFDF2,  (b)  LES-FD.  AE 
=  A,V£=40. 


in  Figs.  6  and  7.  It  is  observed  that  NE  does  not  have  a 
significant  influence  on  the  first  moments,  but  does  slightly 
influence  the  second  moments.  In  all  the  cases  considered, 
NE**40  yields  reliable  predictions,  consistent  with  previous 
consistency  and  convergence  assessments  of  the  scalar 
FDF.6,8  All  the  subsequent  simulations  are  conducted  with 
A£  =  A/2  and  NE=40. 

D.  Comparative  assessments  of  the  VFDF 

The  objective  of  this  section  is  to  analyze  some  of  the 
characteristics  of  the  VFDF  via  comparative  assessments 
against  DNS  data.  This  assessment  is  done  via  both  a  priori 
and  a  posteriori  analyses.  In  the  former,  the  DNS  results  are 
used  to  determine  the  range  of  the  empirical  constants  ap¬ 
pearing  in  the  VFDF  sub-closures.  In  the  latter,  the  final  re¬ 
sults  as  predicted  by  the  VFDF  are  directly  compared  with 
those  obtained  by  DNS.  The  procedure  is  similar  to  that  in 
Ref.  20  and  considers  the  3D  temporal  mixing  layer. 

In  addition  to  VFDF,  three  other  LES  are  conducted  with 
(1)  no  SGS  model,  (2)  the  Smagorinsky 16,60  SGS  closure, 
and  (3)  the  dynamic  Smagorinsky17"19  model.  In  the  case 
with  no  model,  the  contribution  of  the  SGS  is  completely 
ignored,  i.e.,  t£(mi*,m/)  =  0.  In  this  case,  the  numerical  errors 
amount  to  an  implied  model.  But  as  indicated  in  Ref.  20  this 


case  is  included  to  provide  a  point  of  reference  for  the  other 
closures.  The  Smagorinsky  model  is16,61 

Tl(  « I .  Uj )  -  j k  Sij  =  -  2  v,Sij , 


1  /  d(ut)L  djujhX 

iJ  2  \  dxj  dxi  )’ 


(36) 


p,=  C„AlS. 

Cv-y/l  0.172=“0.04,  S=  and  A L  is  the  characteristic 

length  of  the  filter.  This  model  considers  the  anisotropic  part 
of  the  SGS  stress  tensor  atj  =  r£(uz- , uj)  -  \k  8^ .  The  isotro¬ 
pic  components  are  absorbed  in  the  pressure  field.  The  dy¬ 
namic  version  of  the  Smagorinsky  model  provides  a  means 
of  approximating  Cv  as  suggested  in  Refs.  17-19.  The  pro¬ 
cedure  for  the  implementation  of  this  model  in  the  3D  tem¬ 
poral  mixing  layer  LES  is  described  by  Vreman;20  thus  it  is 
not  repeated  here.  (See  Refs.  11,  23,  62,  and  63  for  recent 
reviews  on  SGS  closure  strategy.) 

In  addition  to  the  resolved  velocity  field,  the  primary 
integral  statistical  quantities  considered  for  comparative  as¬ 
sessments  are 


1204  Phys.  Fluids,  Vol.  14,  No.  3,  March  2002 


Gicquel  et  al. 


(a) 


LES-FD 

FJG.  2.  Scatter  plots  of  the  filtered  velocity  field  as  obtained  via  VFDF2  vs 
LES-FD.  (a),  (u)L;  (b),  (u)L.  A£  =  A,  Af£=40. 
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Ej  is  the  kinetic  energy  of  the  resolved  field,  ev  represents 
the  viscous  molecular  dissipation  rate  directly  from  the  fil¬ 
tered  field,  Pk  is  the  production  rate  of  the  SGS  kinetic  en- 
ergy  (or  the  rate  of  energy  transfer  from  the  resolved  filtered 
motion  to  the  SGS  motion),  and  Bk  is  the  total 
backscatter.64"66  The  resolved  molecular  dissipation  rate  is 
always  positive  (by  definition),  but  the  production  rate  of  the 
SGS  kinetic  energy  can  be  locally  negative.  This  backscatter 
is  not  represented  in  the  Smagorinsky  model.  The  dynamic 
model  is  potentially  capable  of  accounting  for  it,  but  at  the 
expense  of  causing  numerical  instabilities.  In  the  implemen¬ 
tation  of  the  dynamic  model  used  here,  backscatter  is 
avoided  by  averaging  the  numerator  and  denominator  of  the 
expression  determining  Cv  (Refs.  19  and  20)  over  the  homo¬ 
geneous  directions.  If  negative  values  are  still  present,  they 
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FIG.  3.  Reynolds  averaged  values  of  the  filtered  streamwise  velocity,  (a) 
Cross-stream  variations  at  *  =  7,  (b)  streamwise  variation  at  y  =  0  (center- 
line).  JV£=40. 


are  set  equal  to  zero.20’63  The  “resolved”  components  of  the 
Reynolds-averaged  stress  tensor  are  denoted  by  Rtj  where 
Rij=((ui)L~(ui)iX(uj)L-(uj)L)-  The  “totaT  Reynolds 
stresses  are  denoted  by  rjj  where  rf</-=(nl—  Ui)(Uj  —  Uj). 
These  are  approximated  by  R +  rL( u z- ,  ujj . 20,67,68  In 
DNS,  the  total  stresses  are  evaluated  directly  and  the  results 
indicate  that  /?y+  rL(ui  ,uj)  does  indeed  approximate 
with  a  maximum  error  of  less  than  10%. 

Figure  8  shows  the  distribution  of  the  particle  number 
density  within  the  whole  computational  domain.  Assuring  an 
approximately  uniform  distribution,  the  values  of  the  mo¬ 
ments  within  local  ensembles  are  compared  with  those  of 
filtered  DNS  data.  These  DNS  data  are  transposed  from  the 
original  high  resolution  1933  points  to  the  low  resolution  of 
333  points,  and  then  are  compared  with  LES  results  on  these 
coarse  points. 

The  DNS  data  are  also  used  to  make  a  priori  estimates 
of  the  model  constants.  The  primary  terms  which  require 
closure  are  the  SGS  dissipation  and  the  velocity-pressure 
scrambling  tensors.  The  model  equation  [Eq.  (20)]  involving 
CE  is  in  a  scalar  form.  For  an  estimate  of  Cl  (thus  C0),  we 
consider  the  following  norm  of  the  corresponding  closure 
[Eq.  (27)3 

II  -  n,7  -(Bij-je  <s,7)|| «  Cj  toll  TL(Ui,Uj)  -  \kSij  II ,  (39) 

where  ||  Wi7||  =  \fW~W~j.  To  estimate  the  coefficients,  a  linear 
regression  is  performed  on  all  the  data  points  at  each  com¬ 
putational  time  step.  The  optimized  constants  as  obtained  in 
this  way  are  denoted  by  Ce  and  Cj.  This  procedure  is  also 
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FIG.  4.  Cross-stream  variations  of  the  Reynolds  averaged  values  of  some  of 
the  components  of  the  SGS  stress  tensor  at  jc  =  7  with  7V£=40.  The  LES-FD 
results  are  obtained  with  A£-0.5A,  NE~  40. 


followed  for  the  Reynolds  averaged  data,  with  the  optimized 
models  obtained  in  this  way  denoted  by  C£  and  Cx.  The 
temporal  variations  of  these  estimated  values  are  shown  in 
Fig.  9.  The  nonuniformity  of  the  coefficients  indicates  the 
“nonuniversality”  of  the  models.  This  is  expected  as  the  flow 
evolves  from  an  initially  smooth  laminar  state  to  a  strong 
three-dimensional  state  (at  /^40)  before  the  action  of  the 
small  scales  becomes  significant.  The  closures  as  adopted  are 
not  fully  suitable  for  application  in  all  of  these  flow  regions. 
Nevertheless,  Fig.  9  indicates  that  the  values  for  these  coef¬ 
ficients  as  suggested  in  RAS,  i.e.,  CX^AA5,  C£^\  are  rea¬ 
sonable,  at  least  within  the  turbulent  regime.  The  influences 
of  these  parameters  are  further  investigated  via  a  posteriori 
analysis  of  the  results  as  discussed  below. 

Figures  10  and  11  show  the  contours  of  the  spanwise  and 
the  streamwise  components  of  the  vorticity  field,  respec¬ 
tively,  at  time  t  —  80.  By  this  time,  the  flow  has  gone  through 
several  pairings  and  exhibits  strong  3D  effects.  This  is  evi¬ 
dent  by  the  formation  of  large  scale  spanwise  rollers  with 
presence  of  counter-rotating  streamwise  vortex  pairs  in  all 
the  simulations.  The  results  via  the  no-model  indicate  too 
many  small-scale  structures  which  clearly  are  not  captured 
accurately  on  the  coarse  grid.  The  amount  of  SGS  diffusion 
with  the  Smagorinsky  model  is  very  significant  at  initial 
times.  Due  to  this  dissipative  characteristics  of  the  model,  the 


FIG.  5.  Percentage  of  the  relative  difference  between  the  L2  norms  of  the 
stresses  as  a  function  of  A£/A.  (a)  x=2.8,  (b)  jc=7,  (c)  x=  11.2. 

predicted  results  are  too  smooth  and  only  contain  the  large 
scale  structures.  The  vortical  structures  as  depicted  by  the 
dynamic  Smagorinsky  and  the  VFDF  are  very  similar  and 
predict  the  DNS  results  better  than  the  other  two  models.  The 
results  obtained  by  VFDF1  and  VTDF2  are  virtually  indis¬ 
tinguishable  from  each  other.  This  is  expected,  due  to  the 
lack  of  importance  of  molecular  effects  in  this  free  shear 
flow. 

The  Reynolds  averaged  values  of  the  streamwise  veloc¬ 
ity  and  the  temporal  variations  of  the  momentum  thickness 

1  CLf 2  -  - 

sm(  0  4  =J_LJ1  “<“>/.)( 1  (4°) 

are  shown  in  Figs.  12  and  13,  respectively.  In  Fig.  12  the 
Reynolds  averaged  values  of  both  filtered  and  unfiltered 
DNS  data  are  considered  and  are  shown  to  be  essentially 
equivalent.  Therefore,  the  latter  are  not  shown  in  subsequent 
figures.  The  dissipative  nature  of  the  Smagorinsky  model  at 
initial  times  resulting  in  a  slow  growth  of  the  layer  is  shown. 
Several  values  of  the  model  parameters  (C0,  C£)  are  consid- 
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FIG.  8.  Particle  number  density  in  VFDF2  simulation  at  r  =  60.  The  isosur¬ 
face  corresponds  to  NE~ 40  set  as  initial  conditions.  C0=2.1,  Ce=  1. 


FIG.  6.  Cross-stream  variations  of  the  Reynolds  averaged  values  of  the 
filtered  streamwise  velocity  at  x-1.  The  LES-FD  results  are  obtained  with 
A£— 0.5A,  Ne=  40. 


y 


y 

FIG.  7.  Cross-stream  variations  of  the  Reynolds  averaged  values  of  some  of 
the  components  of  the  SGS  stress  tensor  at  x-1.  The  LES-FD  results  are 
obtained  with  A£=0.5A,  NE- 40. 


ered  in  the  VFDF  simulations.  It  is  observed  that  as  the  mag¬ 
nitude  of  Ce  decreases,  the  initial  rate  of  the  layer’s  spread  is 
higher.  With  the  exception  of  the  case  with  Ce  =  0.5  and  the 
Smagorinsky  model,  all  the  other  VFDF  cases,  the  dynamic 
Smagorinsky  and  the  no-model  yield  a  similar  rate  of  layer’s 
growth  at  late  times. 

The  temporal  variations  of  the  resolved  kinetic  energy 
and  all  of  the  terms  defined  in  Eq.  (38)  are  shown  in  Fig.  14. 


FIG.  9.  Time  variation  of  the  model  coefficients  as  obtained  from  a  priori 
analysis  of  the  DNS  data. 
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FIG.  10.  Contour  plots  of  the  spanwise  component  of  the  vorticity  at  z 
=  0.7 5LILr,  /=80.  (a)  Filtered  DNS,  (b)  no  model,  (c),  Smagorinsky 
model,  (d)  dynamic  Smagorinsky  model,  (e)  VFDF2,  C0  =  2.1,  C£=l,  (0 
VFDF1,  C0=2.1,  Ce=  1. 


The  overall  features  displayed  in  this  figure  are  similar  to 
those  reported  by  Vreman  et  al20  for  the  no  model,  the  Sma¬ 
gorinsky  model  and  the  dynamic  Smagorinsky  model.  The 
initial  rate  of  decay  of  the  resolved  kinetic  energy  for  the 
Smagorinsky  model  is  the  highest.  This  is  due  to  the  exces¬ 
sive  production  of  the  SGS  kinetic  energy  by  this  model  in 
the  transitional  region,  and  explains  the  reason  for  the  lack  of 
small  scales  in  the  vortical  structures  as  discussed  before.  For 
all  the  other  models  the  initial  rate  of  decrease  of  the  re¬ 
solved  kinetic  energy  is  small  and  increases  as  the  flow  de¬ 
velops.  The  trend  portrayed  by  DNS  results  is  best  captured 
by  the  VFDF  simulations.  For  the  no  model  case  the  only 
means  of  dissipation  of  the  resolved  kinetic  energy  is 
through  molecular  action  and  numerical  dissipation  which 
become  significant  at  later  stages  due  to  presence  of  a  large 
amount  of  small  scales.  In  this  case,  the  amount  of  numerical 
dissipation  is  the  highest.  For  all  the  other  closures,  the  pro¬ 
ductions  rate  of  the  SGS  kinetic  energy  is  larger  than  the 
molecular  dissipation  as  the  flow  develops.  The  dynamic 
Smagorinsky  and  the  no-model  simulations  predict  the  same 
initial  rate  of  decay  for  the  resolved  kinetic  energy.  This  is 
due  to  low  initial  values  of  Pk  predicted  by  the  dynamic 
Smagorinsky  model.  After  f= 40  the  amount  of  Pk  as  pre¬ 
dicted  by  the  dynamic  model  is  more  than  that  of  molecular 
dissipation  by  the  no-model.  Thus  the  rate  of  decay  of  the 
resolved  kinetic  energy  becomes  higher  for  the  dynamic 


(a)  (b) 


FIG.  11.  Contour  plots  of  the  streamwise  component  of  the  vorticity  vector 
at  jc=0.25 L!Lri  /  =  80.  (a)  Filtered  DNS,  (b)  no  model,  (c)  Smagorinsky 
model,  (d)  dynamic  Smagorinsky  model,  (e)  VFDF2,  C0=2.1,  CE=1,  (f) 
VFDF1,  C0  =  2.1,  CE=1. 

model  and  is  closer  to  that  obtained  by  DNS. 

With  the  exception  of  the  no-model  case,  all  the  simula¬ 
tions  predict  similar  trends  for  the  molecular  dissipation.  The 
magnitude  of  this  dissipation  as  predicted  by  VFDF  changes 
slightly  with  the  variation  of  the  model  parameter.  The  pro¬ 
duction  rate  of  the  SGS  kinetic  energy  depends  more 
strongly  on  the  model  coefficients;  as  Ce  decreases,  the  peak 


FIG.  12.  Cross-stream  variations  of  the  Reynolds  averaged  values  of  the 
streamwise  velocity  at  t  =  70. 
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FIG.  13.  Temporal  variations  of  the  momentum  thickness. 


magnitude  of  Pk  is  larger.  The  Smagorinsky  model  does  not 
adequately  predict  Pk ,  and  the  dynamic  model  yields  better 
predictions  at  long  times.  The  overall  trends  are  best  pre¬ 
dicted  by  VFDF.  The  same  is  true  in  capturing  the  backscat- 
ter  phenomenon.  By  design,  the  backscatter  is  identically 
zero  in  the  Smagorinsky  and  the  dynamic  Smagorinsky 
model.  But  VFDF  is  capable  of  capturing  it,  and  its  extent  is 
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o — o 

Dynamic  Smagorinsky 
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■ — ■ 
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FIG.  15.  Cross-stream  variations  of  some  of  the  components  of  atJ  at  t 
=  60. 

controlled  by  the  model  parameters.  In  this  regard  it  is  im¬ 
portant  to  note  that  there  are  no  numerical  instability  prob¬ 
lems  in  the  VFDF  solver  for  negative  Bk  values.  However, 
the  amount  of  predicted  backscatter  is  less  than  that  of 


FIG.  14.  Temporal  variations  of  (a)  total  resolved  ki¬ 
netic  energy,  (b)  SGS  kinetic  energy  production  rate,  (c) 
total  backscatter,  (d)  total  resolved  dissipation. 
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FIG.  16.  Cross-stream  variations  of  some  of  the  components  of  atj  at  t 
=  80. 

DNS  and  its  relative  magnitude  is  less  than  those  of  Pk  and 
Ev. 

Several  components  of  the  planar  averaged  values  of  the 
SGS  anisotropy  tensor,  =  rL{Ui  ,u7)  —  \kStj  are  presented 
in  Figs.  15  and  16.  Both  the  Smagorinsky  and  the  dynamic 
model  under-predict  the  components  of  this  stress.  The 
VFDF  predictions  are  more  satisfactory.  In  this  regard,  the 
VFDF  is  expected  to  be  more  effective  than  the  other  clo¬ 
sures  for  LES  of  reacting  flows  since  the  extent  of  SGS  mix¬ 
ing  is  influenced  by  SGS  convection.69,70  “Optimum”  values 
for  CE  and  C0  cannot  be  suggested  to  predict  all  of  the  com¬ 
ponents  of  this  tensor  at  all  times,  but  it  is  obvious  that  there 
is  too  much  SGS  energy  with  Cfi  =  0.5. 

Several  components  of  the  resolved  stress  tensor  /?/;  are 
shown  in  Figs.  17  and  18.  As  expected,  the  performance  of 
the  Smagorinsky  model  is  not  very  good  as  it  does  not  pre¬ 
dict  the  spread  and  the  peak  value  of  the  resolved  Reynolds 
stresses.  None  of  the  other  models  show  a  distinct  superiority 
in  predicting  the  DNS  results.  The  no-model  and  the  dy¬ 
namic  Smagorinsky  model  predict  large  peak  values  at  the 
middle  of  the  layer.  The  VFDF  predicts  both  the  spread  and 
the  peak  values  reasonably  well.  The  results  for  small  CE 
values  are  not  shown  since  the  amount  of  energy  in  the  re¬ 
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FIG.  17.  Cross-stream  variations  of  some  of  the  components  of  R(j  at  t 
=  60. 

solved  scale  decreases  too  much  in  favor  of  the  increase  of 
the  SGS  stress  (as  shown  in  Figs.  15  and  16).  The  cross¬ 
stream  variations  of  the  total  Reynolds  stress  r12  are  pre¬ 
sented  in  Fig.  19.  The  peak  values  by  the  no-model  simula¬ 
tions  are  again  the  highest.  The  dynamic  model  and  VFDF 
perform  similarly  and  capture  the  DNS  trends  equally  well. 

E.  Comparison  with  previous  investigations 

All  of  the  results  obtained  here  by  DNS,  and  LES  via  the 
Smagorinsky  and  the  dynamic  Smagorinsky  models  agree 
very  well  with  those  of  Vreman  et  al.20  The  slight  differences 
are  due  to  the  nonidentical  flow  initializations,  and  the  dif¬ 
ferent  computational  methodologies  employed  in  the  two 
simulations.  To  compare  with  results  of  other  investigations, 
simulations  are  conducted  of  another  temporally  developing 
mixing  layer  with  Re  =  500  in  a  larger  computational  do¬ 
main,  Lr=  120.  An  initial  forcing  of  the  form  Ae~^y/2)  is 
used,  where  A  is  a  uniformly  distributed  random  number 
with  an  amplitude  of  0.05.  Rogers  and  Moser60  perform  DNS 
of  a  high  Re  number  flow  on512X210Xl92  spectral  points. 
The  results  of  these  simulations  are  in  excellent  agreements 
with  laboratory  data  of  Bell  and  Mehta.71  Here,  LES  is  con 
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FIG.  18.  Cross-stream  variations  of  some  of  the  components  of  /?J;  at  t 
=  80. 

ducted  of  this  flow  via  the  dynamic  Smagorinsky  model. 

The  profiles  of  the  mean  streamwise  velocity  and  several 
components  of  the  resolved  stresses  at  f=250  are  presented 
in  Figs.  20  and  21,  respectively.  In  these  figures,  £ 
—yl <5m(f)  and  the  symbols  denote  the  experimental  data71  at 
several  streamwise  locations.  The  good  agreement  with  these 
data  also  indicates  good  agreement  with  DNS  results  of  Rog¬ 
ers  and  Moser.72 

F.  Computational  requirements 

The  total  computational  times  associated  with  simula¬ 
tions  of  the  3D  temporal  mixing  layer  are  shown  in  Table  EL 
Expectedly,  the  overhead  associated  with  the  VFDF  simula¬ 
tion  is  extensive  as  compared  to  the  other  models;  neverthe¬ 
less  this  requirement  is  significantly  less  that  of  DNS.  This 
overhead  was  tolerated  in  present  simulations,  but  can  be 
reduced  with  utilization  of  an  optimum  parallel  simulation 
procedure.  This  has  been  discussed  for  use  in  PDF73  and  is 
recommended  for  future  VFDF  simulations. 

VII.  SUMMARY  AND  CONCLUDING  REMARKS 

The  filtered  density  function  (FDF)  methodology1  has 
proven  very  effective  for  large  eddy  simulation  (LES)  of 
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FIG.  19.  Cross-stream  variations  of  r12,  (a)  f=60,  (b)  r=80. 

turbulent  reacting  flows.3,6"11  In  all  previous  contributions, 
the  LES/FDF  of  only  the  scalar  quantities  are  considered. 
The  objective  of  the  present  work  is  to  develop  the  FDF 
methodology  for  LES  of  the  velocity  field.  For  this  purpose, 
a  methodology  termed  the  velocity  filtered  density  function 


FIG.  20.  Cross-stream  variation  of  the  Reynolds  averaged  values  of  the 
streamwise  velocity  at  /  =  250.  Solid  line  denotes  model  predictions  via  the 
dynamic  Smagorinsky  model.  Symbols  denote  experimental  data  of  Bell  and 
Mehta  (Ref.  71). 
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FIG.  21.  Cross-stream  variations  of  the  Reynolds  aver¬ 
aged  values  of  the  streamwise  velocity  at  t  —  250.  Solid 
lines  denote  model  predictions  via  the  dynamic  Smago- 
rinsky  model.  Symbols  denote  experimental  data  of 
Bell  and  Mehta  (Ref.  71). 


(VFDF)  is  developed.  The  VFDF  is  basically  the  probability 
function  (PDF)  of  the  subgrid  scale  (SGS)  velocity  vector. 
The  exact  transport  equation  governing  the  evolution  of  the 
VFDF  is  derived.  It  is  shown  that  the  effects  of  SGS  convec¬ 
tion  in  this  equation  appears  in  a  closed  form.  The  unclosed 
terms  in  this  transport  equation  are  modeled  via  two  formu¬ 
lations:  VFDF1  and  VFDF2.  The  primary  difference  between 
the  two  models  is  the  inclusion  of  the  molecular  diffusion  in 
the  spatial  transport  of  the  VFDF  in  the  first  formulation.  The 
closure  strategy  in  the  formulation  similar  to  that  in  PDF 
methods  in  Reynolds  averaged  simulation  (RAS) 
procedures.32  In  this  way,  the  VFDF  formulation  is  at  least 
equivalent  to  a  second-order  moment  SGS  closure. 

The  modeled  VFDF  transport  equations  are  solved  nu¬ 
merically  via  a  Lagrangian  Monte  Carlo  scheme  in  which  the 
solutions  of  the  equivalent  stochastic  differential  equations 
(SDEs)  are  obtained.  Two  Monte  Carlo  procedures  are  con¬ 
sidered.  The  schemes  preserve  the  Ito-Gikhman  nature  of 
the  SDEs  and  provide  a  reliable  solution  for  the  VFDF.  The 
consistency  of  the  VFDF  formulation  and  the  convergence  of 
its  Monte  Carlo  solutions  are  assessed.  This  is  done  via  com¬ 
parisons  between  the  results  obtained  by  the  Monte  Carlo 
procedure  and  the  finite  difference  solution  of  the  transport 
equations  of  the  first  two  filtered  moments  of  VFDF  (LES- 
FD).  With  inclusion  of  the  third  moments  from  the  VFDF 
into  the  LES-FD,  the  consistency  and  convergence  of  the 


TABLE  II.  Computer  requirements  for  the  3D  temporal  mixing  layer.  One 
unit  corresponds  to  1657.2  seconds  of  CPU  time  on  the  SGI  origin  2000. 


Resolution 

Ne 

Normalized  CPU  time 

DNS 

193X193X193 

178 

VFDF1 

33X33X33 

40 

33.6 

VFDF2 

33X33X33 

40 

30 

Dynamic  Smagorinsky 

33X33X33 

2.19 

Smagorinsky 

33X33X33 

1.05 

No  model 

33X33X33 

1 

Monte  Carlo  solution  is  demonstrated  by  good  agreements  of 
the  first  two  SGS  moments  with  those  obtained  by  LES-FD. 

The  VFDF  predictions  are  compared  with  those  with 
LES  results  with  no  SGS  model,  with  the  Smagorinsky16 
SGS  closure,  and  with  the  dynamic  Smagorinsky17-19  model. 
All  of  these  results  are  also  compared  with  direct  numerical 
simulation  (DNS)  results  of  a  three-dimensional,  temporally 
developing  mixing  layer  in  a  context  similar  to  that  con¬ 
ducted  by  Vreman  et  al20  This  comparison  provides  a  means 
of  examining  some  of  the  trends  and  overall  characteristics 
as  predicted  by  LES.  It  is  shown  that  the  VFDF  performs 
well  in  predicting  some  of  the  phenomena  pertaining  to  the 
SGS  transport.  The  magnitude  of  the  SGS  Reynolds  stresses 
as  predicted  by  VFDF  is  larger  than  those  predicted  by  the 
other  SGS  models  and  much  closer  to  the  filtered  DNS  re¬ 
sults.  The  temporal  evolution  of  the  production  rate  of  the 
SGS  kinetic  energy  is  predicted  well  by  VFDF  as  compared 
with  those  via  the  other  closures.  The  VFDF  is  also  capable 
of  accounting  the  SGS  backscatter  without  any  numerical 
instability  problems,  although  the  level  predicted  is  substan¬ 
tially  less  than  that  observed  in  DNS. 

The  results  of  a  priori  assessment  against  DNS  data  in¬ 
dicates  that  the  values  of  the  model  coefficients  as  employed 
in  VFDF  (C0  and  Ce)  are  of  the  range  suggested  in  the 
equivalent  models  previously  used  in  RAS.  The  results  of  a 
posteriori  assessments  via  comparison  with  DNS  data  does 
not  give  any  compelling  reasons  to  use  values  other  than 
those  suggested  in  RAS,  C0  =  2.1,  Ce=l.  However,  small 
values  of  Ce  are  not  acceptable  as  they  would  yield  too  much 
of  SGS  energy  relative  to  that  within  the  resolved  scales. 

Most  of  the  overall  flow  features,  including  the  mean 
velocity  field  and  the  resolved  and  total  Reynolds  stresses  as 
predicted  by  VFDF  are  similar  to  those  obtained  via  the  dy¬ 
namic  Smagorinsky  model.  This  is  interesting  in  view  of  the 
fact  that  the  model  coefficients  in  VFDF  are  kept  fixed.  It 
may  be  possible  to  improve  the  predictive  capabilities  of  the 
VFDF  by  two  ways:  (1)  Development  of  a  dynamic  proce- 
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dure  to  determine  the  model  coefficients,  and/or  (2)  imple¬ 
mentation  of  higher  order  closures  for  the  generalized  Lange- 
vin  model  parameter  Giy  (see  Ref.  34). 

Work  is  in  progress  towards  developments  of  a  joint 
velocity-scalar  FDF  for  LES  of  reacting  flows.  Compared  to 
standard  LES,  this  approach  has  the  advantage  of  treating 
reaction  in  a  closed  form;  and,  compared  to  scalar  FDF6,8  has 
the  advantage  of  treating  convective  transport  (of  momentum 
and  species)  in  closed  form.  These  modeling  advantages 
have  an  associated  computational  penalty.  For  the  cases  con¬ 
sidered  here,  VFDF  is  more  expensive  computationally  than 
the  dynamic  Smagorinsky  model  by  a  factor  of  15.  It  is  ex¬ 
pected  that  VFDF  will  not  be  more  expensive  than  scalar 
FDF,  at  least  for  reacting  flows  with  many  species. 
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A  stochastic  model  is  developed  for  the  acceleration  of  a  fluid  particle  in  anisotropic  and 
inhomogeneous  turbulent  flows.  The  model  consists  of  an  ordinary  differential  equation  for  velocity 
(which  contains  directly  the  acceleration  due  to  the  mean  and  rapid  pressure  gradients),  and  a 
stochastic  model  for  the  remainder  of  the  acceleration,  which  is  due  to  the  slow  pressure  gradient 
and  to  viscosity.  In  addition  to  a  rapid-pressure  model,  the  stochastic  model  involves  three  tensor 
coefficients.  For  isotropic  turbulence,  the  model  reverts  to  that  previously  proposed  by  Sawford.  At 
high  Reynolds  number  the  model  is  consistent  with  local  isotropy  and  the  Kolmogorov  hypotheses, 
and  tends  to  the  generalized  Langevin  model  for  fluid-particle  velocity.  In  this  case  two  of  the  tensor 
coefficients  are  known  in  terms  of  the  Kolmogorov  constant  C0 ,  while  the  third  is  related  to  the 
coefficient  in  the  generalized  Langevin  model.  A  complete  analysis  of  the  model  is  performed  for 
homogeneous  turbulent  shear  flow,  for  which  there  are  Lagrangian  data  from  direct  numerical 
simulations.  The  main  result  is  to  establish  the  one-to-one  correspondence  between  the  model 
coefficients  and  the  primary  statistics,  namely,  the  velocity  and  acceleration  covariances  and  the 
tensor  of  velocity  integral  time  scales.  The  autocovariances  of  velocity  and  acceleration  obtained 
from  the  model  are  in  excellent  agreement  with  the  direct  numerical  simulation  (DNS)  data.  Future 
DNS  studies  of  homogeneous  turbulence  can  be  used  to  investigate  the  dependence  of  the  model 
coefficients  on  Reynolds  number  and  on  the  imposed  mean  velocity  gradients.  The  acceleration 
model  can  be  used  to  generate  a  range  of  turbulence  models  which,  in  a  natural  way,  incorporate 
Reynolds-number  effects.  ©  2002  American  Institute  of  Physics.  [DOI:  10.1063/1.1483876] 


I.  INTRODUCTION 

In  order  to  investigate  dispersion  in  turbulent  flows,  in 
1921  Taylor1  introduced  a  stochastic  model  for  the  position 
X+(r)  of  a  fluid  particle.  An  analysis  of  Taylor’s  model 
shows  that  it  is  equivalent  to  the  Langevin  equation  as  a 
model  for  the  fluid-particle  velocity  U+(f)  =  dX+(t)/dt. 
(Langevin2  had  proposed  this  stochastic  equation  in  1908  to 
model  the  velocity  of  particles  undergoing  Brownian  mo¬ 
tion.)  The  Langevin  equation  remains  the  basis  for  stochastic 
models  of  turbulent  dispersion  (see,  e.g..  Refs.  3-5).  Further¬ 
more  the  Langevin  equation  and  its  generalization6,7  provide 
a  closure  to  the  transport  equation  for  the  (one-point,  one¬ 
time)  probability  density  function  (PDF)  of  velocity.8,9  And 
from  the  modeled  velocity  PDF  equation  can  be  deduced  the 
corresponding  partially  modeled  Reynolds-stress  equation.10 
Thus,  an  accurate  stochastic  model  for  the  fluid-particle 
velocity  U+(/)  is  a  potent  tool  in  turbulence  modeling  as 
well  as  in  the  study  of  turbulent  dispersion. 

Important  conclusions  about  the  performance  of  the 
Langevin  model  can  be  drawn  from  the  simplest  case  of 
statistically  stationary  homogeneous  isotropic  turbulence.  In 
general,  the  fluctuating  component  of  fluid-particle  velocity 
is  defined  by 

u+(r)  =  U+(/)-(U(X+[r],r)),  (1) 

where  U(x,/)  is  the  Eulerian  velocity;  and  for  the  case  con- 
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sidered  u+(r)  is  a  statistically  stationary  process  with  mean 
zero.  The  Lagrangian  velocity  autocorrelation  function  is  de¬ 
fined  by 

P(i)s(«(/)(f)«(i)((+j))/<«j)(')«(;)(f)>,  (2) 

which  is  independent  of  t  and  i  because  of  stationary  and 
isotropy,  respectively.  (Here  and  below,  bracketed  suffixes 
are  excluded  from  the  summation  convention.)  The  Langevin 
model  predicts  this  autocorrelation  function  to  be9 


\  lL  f 

where  TL  is  the  Lagrangian  integral  time  scale.  For  not-too- 
small  time  intervals  \s\/TL,  this  prediction  is  in  excellent 
agreement  with  experimental  and  direct  numerical  simula¬ 
tion  (DNS)  data.11 

But  the  form  of  Eq.  (3)  reveals  three  related  shortcom¬ 
ings  of  the  Langevin  model.  First,  it  contains  the  single  time 
scale  Tl  (which  is  characteristic  of  the  large-scale,  energy- 
containing  motions);  second,  there  is  no  dependence  on  Rey¬ 
nolds  number,  and,  third,  the  slope  of  p(s)  given  by  Eq.  (3) 
is  discontinuous  at  the  origin  [reflecting  the  fact  that  the 
Langevin  model  for  u+(0  is  continuous  but  not  differen¬ 
tiable].  The  same  observations  can  be  made12  regarding  the 
Lagrangian  velocity  frequency  spectrum  EL{  co) — which  is 
the  Fourier  transform  of  (u^u^) p(s) .  According  to  the 
Langevin  model,  at  high  frequency  El(<d)  varies  at  to"2: 
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there  is  no  representation  of  the  more  rapid  decrease  in 
El(oj)  beyond  the  frequency  corresponding  to  the  Kolmog¬ 
orov  time  scale  tv. 

In  1991,  Sawford12  introduced  (for  isotropic  turbulence) 
a  stochastic  model  for  the  fluid-particle  acceleration  A+(/) 
=  dV+(t)/dt-d2X+(t)ldt2.  Such  a  model  remedies  the 
above-mentioned  deficiencies  of  the  Langevin  model:  a  sec¬ 
ond  time  scale  (which  scales  with  rv)  is  introduced;  there  is 
an  intrinsic  Reynolds-number  dependence  (since  TLfrv  in¬ 
creases  with  Reynolds  number);  and,  at  the  origin,  the  pre¬ 
dicted  velocity  autocorrelation  function  is  once  continuously 
differentiable.  Correspondingly,  around  the  Kolmogorov  fre¬ 
quency  r“!,  the  Lagrangian  velocity  spectrum  El(cd) 
smoothly  changes  its  power-law  behavior  from  co~2  to  co~4. 
For  isotropic  turbulence,  Sawford's  model  is  in  excellent 
agreement  with  DNS  data,  including  accounting  for  the 
Reynolds-number  dependence  of  the  acceleration  autocorre¬ 
lation  function  and  the  second-order  Lagrangian  velocity 
stmcture  function.  ’ 

In  this  paper  we  consider  a  more  general  stochastic 
model  for  the  fluid-particle  acceleration,  which  is  applicable 
to  anisotropic  turbulence  and  to  inhomogeneous  turbulent 
flows.  The  general  form  of  the  model  is  developed  in  Sec.  II, 
where  particular  attention  is  paid  to  the  contribution  to  ac¬ 
celeration  from  the  rapid  pressure  gradient.  When  applied  to 
homogeneous  turbulence  (with  constant  and  uniform  mean 
velocity  gradients)  the  stochastic  model  is  of  the  form 

where  u*(/)  is  the  model  for  u+(0,  a*(f)  is  its  rate  of 
change  (i.e.,  a*^^*/^),  and  W (t)  is  an  isotropic  Wiener 
process.9  The  coefficients  B,  C,  and  D  are  tensors  which  can 
depend  on  the  local  state  of  the  flow  and  the  turbulence,  but 
are  independent  of  a*  and  u*.  (The  conventional  notation  is 
that  “  +  ”  denotes  a  fluid-particle  property,  and  denotes  a 
model  for  that  property.) 

In  the  simplest  case  of  isotropic  turbulence,  all  the  coef¬ 
ficients  in  Eq.  (4)  are  isotropic  (e.g.,  B{j  =  B  <S/;),  and  the 
model  reverts  to  Sawford's.12  In  this  case,  which  is  reviewed 
in  Sec.  Ill  C,  there  is  a  one-to-one  correspondence  between 
the  three  scalar  coefficients  ( B ,  C,  and  D)  and  the  three  pri¬ 
mary  statistics:  the  acceleration  variance  a'2;  the  velocity 
variance  w'2;  and  the  velocity  integral  time  scale  TL . 

Beyond  isotropic  turbulence,  the  simplest  type  of  flow  to 
study  is  statistically  stationaiy  homogeneous  turbulence  with 
imposed  mean  velocity  gradients — as  exemplified  by  a  re¬ 
cent  DNS  of  forced  homogeneous  turbulent  shear  flow,14  and 
described  in  Sec.  IV.  For  this  case  the  coefficients  B,  C,  and 
D  are  constant,  and  a  complete  analysis  of  the  model  [Eq. 
(4)]  can  be  performed.  This  is  done  in  Sec.  IV  B,  where  it  is 
shown  that  there  is  a  one-to-one  correspondence  between  the 
tensor  coefficients  in  the  model  and  the  primary  statistics, 
namely  the  velocity-acceleration  covariances  and  the  veloc¬ 
ity  integral  time  scale  tensor.  (This  analysis  parallels  the  au¬ 
thors's  recent  analysis  of  a  stochastic  model  for  velocity.15) 

With  some  approximation  (and  with  an  appropriate  scal¬ 
ing  of  the  variables),  the  same  analysis  can  be  applied  to 
nonstationary  homogeneous  turbulence,  in  particular  to  (un¬ 
forced)  homogeneous  turbulent  shear  flow  for  which  there 


are  Lagrangian  data  from  the  recent  DNS  studies  of  Sawford 
and  Yeung.16,17  It  is  shown  (in  Sec.  IV  C)  that  the  velocity- 
acceleration  autocorrelation  functions  predicted  by  the  model 
are  in  excellent  agreement  with  these  DNS  data. 

As  well  as  being  useful  in  its  own  right,  we  also  regard 
the  acceleration  model  as  an  intermediate  step  in  the  devel¬ 
opment  of  improved  stochastic  models  for  velocity  for  use  in 
dispersion  studies,  in  PDF  methods,  and  in  other  turbulence 
models.  Compared  to  the  velocity  model,  the  acceleration 
model  can  be  more  closely  related  to  Lagrangian  data  from 
DNS,  which  are  known  to  contain  strong  Reynolds-number 
dependencies.11,18  Given  an  acceleration  model  (i.e.,  a  pre¬ 
scription  for  the  coefficients  B,  C,  and  D),  a  corresponding 
velocity  model  can  be  deduced15  which  yields  the  same  ve¬ 
locity  covariances  and  integral  time  scales,  and  which  inher¬ 
its  Reynolds-number  dependencies.  Such  an  improved  model 
has  direct  application  in  PDF  methods,  and  from  it  can  be 
deduced  a  pressure-rate-of-strain  model  for  use  in  Reynolds- 
stress  models.  These  and  other  uses  of  the  acceleration  model 
are  discussed  in  Sec.  V. 


II.  STOCHASTIC  MODEL  FOR  ACCELERATION 


We  consider  the  inhomogeneous  turbulent  flow  of  a 
constant-property  Newtonian  fluid  (of  density  p  and  kine¬ 
matic  viscosity  v).  This  is  governed  by  the  continuity  equa¬ 
tion  dU i/dXi  —  0,  and  the  Navier-Stokes  equation 


Aj(xj)  = 


DU  t 
~dT 


Ui 


1  dp  d2Ui 

p  dX  i  dXjdXj  ’ 


(5) 


where  A(x,/),  U(x,/),  and  p{x,t)  are  the  acceleration,  ve¬ 
locity,  and  pressure.  The  general  fluid  particle  has  position 
X+(f),  velocity, 


U +(t) 


dX+(t) 

dt 


u(X+[?],/), 


(6) 


and  acceleration 

nfU+(t) 

A+(/)=  =A(X  [t],r). 


(7) 


A.  Decomposition  of  acceleration 

The  acceleration  can  be  decomposed  into  mean  and  fluc¬ 
tuating  contributions  based  on  the  mean  «U)  and  (/?))  and 
fluctuating  (u  and  p')  components  of  velocity  and  pressure. 
Furthermore,  as  originally  shown  by  Chou,19  the  fluctuating 
pressure  can  be  decomposed  into  rapid,  p{r\  slow,  pis\  and 
harmonic,  p^h\  contributions.9  Thus,  the  fluid  acceleration  is 

1  d(p)  1  dp(r)  1  dp^_  1  dp(h) 

1  p  dx j  p  dX;  p  dXj  p  dXj 

*u, 

V  dXjdXj  dXjdXj * 

The  harmonic  pressure  and  the  mean  viscous  term  are  neg¬ 
ligible  except  in  the  immediate  vicinity  of  walls  (or  other 
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surfaces).  Here  we  neglect  these  terms,  and  hence  leave  to 
future  work  the  development  of  the  special  treatments  re¬ 
quired  for  the  viscous  near-wall  region. 


B.  Structure  of  the  model 


The  proposed  model  consists  of  an  ordinary  differential 
equation  (ODE)  for  U*(t) — a  model  for  the  fluid  particle 
velocity  U+(/) — and  a  stochastic  differential  equation  (SDE) 
for  an  acceleration  variable  denoted  by  A°(r).  The  model 
also  involves  the  fluctuating  components  of  these  quantities, 
which  are  defined  by 


u*(/)-U*(/)-<U*(/)|X*(f)> 

(9) 

a0(0— A°(/)— (A°(f)|X*(/)>, 

(10) 

where  X*(t)  denotes  the  position  of  the  model  particle. 
The  ODE  for  velocity  is 


dUf 

~dT 


\  <?(/>)  \ 
p  dXj  j  x*(,) 


1 

P  J 


+  «?(*), 

x*(/) 


(11) 


where  the  first  two  contributions  on  the  right-hand  side  rep¬ 
resent  acceleration  by  the  mean  pressure  gradient  (which  is 
assumed  to  be  known),  and  acceleration  by  the  rapid  pressure 
gradient  (which  has  to  be  modeled).  A  comparison  of  Eq.  (8) 
and  Eq.  (11)  then  reveals  that  a°(t)  is  a  model  for  the  accel¬ 
eration  due  to  the  slow  pressure  gradient  and  the  viscous 
term. 

The  acceleration  variable  A°(r)  is  modeled  by  the  gen¬ 
eral  SDE, 


dA<?(t)=-[CijA?(t)  +  Dijuf(t)]dt  +  BijdWj,  (12) 

where  W(f)  is  an  isotropic  Wiener  process.  The  tensor  func¬ 
tions  B(x,0,  C(x,/),  and  D(x,0  [which  in  Eq.  (12)  are 
evaluated  at  [X*(/L),t]]  depend  on  the  local  state  of  the  tur¬ 
bulence,  but  are  independent  of  U*  and  A0. 


C.  Homogeneous  turbulence 

Before  presenting  the  rationale  for  the  structure  of  the 
model,  we  first  note  the  form  that  it  takes  in  homogeneous 
turbulence. 

In  homogeneous  turbulence  (with  uniform  mean  velocity 
gradients),  the  coefficients  B,  C,  and  D  depend  only  on  time, 
and  it  follows  that  the  mean  (A0(f)|X*(O)  =  (A0(/))  is  zero. 
Consequently  a°(/)  is  identical  to  A°(r).  And  the  velocity 
equation  can  readily  be  transformed  to  an  equation  for 
u*^)-  Thus,  for  homogeneous  turbulence  the  model  be¬ 
comes 


sure  gradients — appear  directly  in  the  ODE  for  velocity,  Eq. 
(11);  whereas  the  other  contributions — from  the  slow  pres¬ 
sure  gradient  and  viscosity — are  modeled  through  the  SDE 
for  A°(0,  Eq.  (12).  The  rationale  for  this  division  is  based 
on  the  response  of  the  system  to  a  rapid  distortion,  and  it  can 
be  most  easily  understood  for  the  case  of  homogeneous  tur¬ 
bulence. 

Consider  the  sudden  imposition  of  a  very  large  strain 
rate  on  homogeneous  turbulence.  Both  the  mean  and  rapid 
pressure  fields  change  suddenly  and  this  leads  to  a  sudden 
change  in  the  fluid  acceleration.  On  the  other  hand,  the  fluc¬ 
tuating  velocity  field  and  the  slow  pressure  change  continu¬ 
ously  in  response  to  the  suddenly  imposed  distortion. 

The  model  is  qualitatively  in  accord  with  this  behavior. 
It  may  be  seen  from  Eq.  (11)  and  Eq.  (13)  that  the  accelera¬ 
tion  changes  suddenly  if  there  is  a  sudden  change  in 
d{U)ldxj,  with  accompanying  sudden  changes  in  d(p)fdXi 
and  dp{r)/dXi .  In  the  acceleration  equation  [Eq.  (12)  and  Eq. 
(14)],  these  sudden  changes  can  result  in  sudden  changes  in 
the  coefficients,  B,  C,  and  D,  but  nevertheless,  a °(f)  changes 
continuously. 


E.  Rapid-pressure  models 


As  is  usual,  and  in  keeping  with  the  physics,  we  consider 
deterministic  models  for  the  rapid  pressure  gradient.  The 
quantity  then  to  be  modeled  is  the  conditional  mean  rapid 
pressure  gradient — conditional  on  the  modeled  state  of  the 
fluid  particle. 

For  homogeneous  turbulence,  the  rapid  pressure  varies 
linearly  with  the  imposed  mean  velocity  gradient,  and  hence 
the  general  model  can  be  written20 


1  dp(r) 
p  dx{ 


=  2 


dxt 


AW 


(15) 


The  third-order  tensor  function  Neki  is  given  in  terms  of 
two-point  conditional  velocity  statistics  in  Refs.  20  and  8 
(where  it  is  denoted  by  Beki),  and  it  satisfies  the  relations 


A reu—u*’  N€€i— 0,  N€ki—Neik.  (16) 

Rapid  distortions  of  homogeneous  turbulence  can  be 
treated  exactly  using  the  wave-vector  model  of  Van  Slooten 
and  Pope.9,21  This  requires  that  the  modeled  state  of  the  fluid 
particle  be  supplemented  by  the  wave  vector  e*(f) — which, 
among  other  conditions,  satisfied  the  relations 

ef  ef=  1,  -  ef  uf  —  0.  (17) 

Then  the  tensor  Neki  in  Eq.  (15)  is  given  by 

N(ki=ueeke?  ■  (18) 


duf 

~dT 


—z - uf- 

dx :  J 


*  *P(r) \ 
p  dxi  j 


+«?(*). 


x*(0 


da%t)  =  —[C jja'jit)  +  D ijUj  (t)]dt + B ij  dWj . 


(13) 

(14) 


D.  Rationale 

The  structure  of  the  model  is  such  that  some  contribu¬ 
tions  to  acceleration — namely,  from  the  mean  and  rapid  pres- 


For  rapid  distortions,  the  wave-vector  model  consists  of 
ODE's  for  e*(t)  and  u*(f),  the  latter  being  Eq.  (13)  with  the 
neglect  of  a0,  and  with  the  rapid-pressure  model  given  by 
Eqs.  (15)  and  (18).  This  model  is  exact  for  arbitrary  rapid 
distortions  of  homogeneous  turbulence,  in  the  sense  that  it 
yields  the  correct  evolution  of  the  Reynolds  stresses. 

As  is  conventional  in  Reynolds-stress  and  velocity-PDF 
modeling,  we  are  primarily  concerned  here  with  models 
based  on  velocity  and  its  one-point  statistics,  i.e.,  u*(f)  and 
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the  Reynolds  stress  (w/wy).  The  unfortunate  fact  of  the  matter 
is  that  these  quantities  are  inadequate  to  describe  rapid  dis¬ 
tortions  (see,  e.g.,  Reynolds  and  Kassinos22):  additional  di¬ 
rectional  information  is  needed,  as  is  provided  by  the  wave 
vector.  However,  the  hope  is  that  rapid-pressure  models 
based  on  velocity  alone  may  be  adequate  for  the  moderate 
and  slowly  varying  mean  strain  rates  encountered  in  many 
turbulent  shear  flows. 

Following  Ref.  20,  it  is  natural  to  consider  a  rapid- 
pressure  model  that  is  linear  in  velocity,  and  which  therefore 
can  be  written 


i  dp{r) 
P  dxt 


0  jfc 

a  ,u* 


=//<>; 


dxe  ’ 


(19) 


III.  PROPERTIES  OF  THE  MODEL 

In  this  section  we  examine  some  of  the  mathematical 
properties  of  the  model,  and  their  connections  to  the  physics 
of  turbulent  motions. 

A.  Equivalent  first-order  and  second-order  systems 

For  homogeneous  turbulence,  the  model  [Eq.  (14)  and 
Eq.  (22)]  can  be  written  as  a  first-order  system  of  SDE’s, 


duf=[-Kyuf+a^]dt, 

(25) 

da<>=  -[ Cya^+DyufJdt+BydWj , 

(26) 

or,  in  an  inferior  notation,  as  a  first-order  system  of  ODE’s 


or,  equivalently, 

N  ai~  lH\jleuf »  @0) 

where  the  tensors  G-J*  and  Hyle  correspond  to  the  analogous 
tensors  in  the  Haworth-Pope  model.7,9  The  nondimensional 
fourth-order  tensor  H(r>  is  modeled  as  a  linear  function  of  the 
Reynolds-stress  anisotropy  tensor 


(UjUj) 
ij  (Uk“k) 


(21) 


and  indeed  a  nontrivial  dependence  on  by  is  required  to  sat¬ 
isfy  the  condition  that  the  rapid  pressure  neither  produces  nor 
removes  turbulent  kinetic  energy. 

In  subsequent  sections  we  confine  attention  to  this  linear 
model,  not  least  because  it  is  amenable  to  analysis.  In  DNS, 
the  rapid  pressure  gradient  can  extracted,  its  linearity  in  u+ 
can  be  examined,  and  specific  models  for  Hyl(  can  be  as¬ 
sessed.  With  this  model,  Eq.  (13)  can  be  rewritten 


duf 

~dt~ 


=  -KiJuf+a", 


(22) 


d±  =  -K, juf+J, 


(27) 


da 


dt 


L  —  Ctri-Dijuf  +  ByWj, 


(28) 


where  W  denotes  white  noise,  which  has  the  property 

JoW(f')df'  =  W(f). 

Alternatively,  by  differentiating  Eq.  (27)  with  respect  to 
t,  the  model  can  be  re-expressed  as  the  second-order  system 


2„* 


dzu 


du' 


dt 


r  +  (Cij+Kij)-jf-+  Djj+ 


dKj 

~dt 


uf  =  BijWj.  (29) 


It  may  be  seen  that  the  system  is  governed  fundamentally  by 
just  three  coefficient  tensors,  not  four  as  suggested  by  the 
appearance  of  B,  C,  D,  and  K  in  Eqs.  (25)-(28).  In  particu¬ 
lar,  if  K  is  constant — as  is  the  case  in  the  analysis  below 
(Sec.  IV  B) — the  behavior  of  u*(r)  is  determined  by  B,  D 
and  the  sum 


C=C+K, 


(30) 


where  the  tensor  Ky  is  defined  by 


Ky  = 


dxj  dX(> 


Wjt-rflti 


(23) 


Finally,  we  caution  that  in  future  studies  of  rapid- 
pressure  modeling  (e.g.,  based  on  DNS)  nonlinear  models 
should  not  be  discounted.  For  example,  the  simple  model 


%=2Wf|  Si> 


ik  ufuf 


(24) 


satisfies  all  known  constraints  (without  requiring  a  depen¬ 
dence  on  bjj). 


F.  Summary 

The  model  consists  of  an  ODE  for  velocity,  Eq.  (11), 
which  contains  the  mean  pressure  gradient  and  a  model  for 
the  rapid  pressure  gradient  [e.g.,  Eq.  (19)].  The  remainder  of 
the  acceleration — owing  to  the  slow  pressure  gradient  and 
the  viscous  term — is  modeled  by  an  SDE,  Eq.  (12),  which 
contains  three  tensor  coefficients,  B,  C,  and  D.  Various 
properties  of  these  coefficients  are  revealed  in  subsequent 
sections. 


but  not  by  C  and  K  individually.  Thus,  for  constant  K,  Eqs. 
(25)  and  (26)  are  equivalent  to  the  system 

duf  =  afdt ,  (31) 

daf  =  ~[Cuaf  +  Duuf]dt + Bu  dWj ,  (32) 

in  which  Eq.  (31)  defines  n*=du*ldt,  and  a0  can  be  recov¬ 
ered  as 

a^af  +  Kijuf.  (33) 

The  model  is  analyzed  below  via  Eqs.  (31)  and  (32). 

B.  Scaled  model  equations  and  coefficients 

It  is  informative  to  scale  the  variables  and  coefficients  in 
the  model  equations  for  homogeneous  turbulence  so  that 
they  become  nondimensional  quantities  of  order  unity.  A  pre¬ 
liminary  is  to  define  the  quantities  used  to  perform  these 
scalings. 

The  velocity  and  acceleration  variables  are  scaled  by 
their  standard  deviations  uf  and  a\  which  are  given  by 

a’^afaf),  (34) 
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TABLE  I.  Summary  of  different  time  scales. 


T^kfe 

Turbulence  time  scale 

Tv^{vle)m 

Kolmogorov  time  scale 

TS=S~{ 

Shear  time  scale 

in 

Lagrangian  velocity  integral  time  scale 

Ta^u,2!(a'2T) 

Acceleration  time  scale 

Eq.  (55) 

Velocity  eigen-time  scale 

To — ^  7 1 .  Eq.  (57) 

Acceleration  eigen-time  scale 

T,  Eq.  (102) 

Integral  time  scale  (6x6)  matrix 

Tu“,  Eq.  (109) 

Velocity  integral  time  scale  tensor 

where  k  is  the  turbulent  kinetic  energy.  There  are  four  rel¬ 
evant  time  scales.  The  turbulence  time  scale  is  defined  by 


The  ODE  for  velocity,  Eq.  (13)  and  Eq.  (15),  can  be 
written 


1  duf  _  t  duf 
u'  dt  u'  dt 


r\  d(Uk) 
—  Tg-Z - 

TsL  9xt 


2Neki 

u' 


u* 


(43) 


Each  term  is  nondimensional;  expressions  in  square  brackets 
are  of  order  unity;  and  time  is  scaled  by  the  turbulence  time 
scale,  i.e.,  t=t/r.  If  the  linear  rapid-pressure  model,  Eq. 
(19),  is  used,  then  the  ODE  for  u*(f)  can  alternatively  be 
written 


(35) 


where  €  is  the  rate  of  dissipation  of  k.  The  shear  time  scale, 
characteristic  of  the  imposed  mean  velocity  gradients,  is  de¬ 
fined  by 

rs=S~\  (36) 

where 


^2_  d(u,)  d(u,) 

dXj  dXj 

The  Kolmogorov  time  scale  is 


and  the  acceleration  time  scale  is  defined  by 


(37) 


(38) 


(39) 


The  ratio  tv!t  decreases  with  Reynolds  number  as 


1/2 


r: 


(40) 


1  duf 

7d~dT 


T  ~ 
~KU 

7s  J\ 


T 

+  1  — 
T. 


1/21 


(44) 


where  the  nondimensional,  order-one  coefficient  K  is 


Ku = rsKir  Ts-^-[SikSj€ - H%],  (45) 

It  is  clear  from  Eqs.  (43)  and  (44)  that  u flu9  responds  to 
the  mean  velocity  gradients  at  the  normalized  rate  t!ts 
~  Ski  6.  Under  usual  circumstances  this  is  of  order  one,  but 
for  rapid  distortions  it  is  arbitrarily  large.  Evidently,  the  term 
in  a0  is  of  order  yjr/ra^  Re1/4.  But  since  a0  is  a  zero-mean 
random  function  with  normalized  time  scale  Tafr,  the  cumu¬ 
lative  effect  of  the  term  on  the  covariances  of  u */u '  over  a 
time  interval  A i>ra/T  is  of  order  \jrl ra2(\ira It) —At. 
Thus,  although  the  term  in  a0  is  relatively  large  instanta¬ 
neously  (of  order  Re1/4),  its  cumulative  effect  is  of  order  one. 

For  the  SDE  for  a°(f),  Eq.  (14),  we  define  the  scaled 
coefficients  by 


C=raC,  b=rraD. 


(46) 


where  the  turbulence  Reynolds  number  is  R t=k?-/(€v),  and 
the  Taylor-scale  Reynolds  number  is  Rk  =  ( yRe)I/2.  The 
various  time  scales  used  throughout  the  paper  are  summa¬ 
rized  in  Table  I. 

With  a0  being  the  Kolmogorov-scaled  acceleration  vari¬ 
ance 

a'2r 

(41) 


the  acceleration  and  Kolmogorov  time  scales  are  related  by 


7a 


T 


V 


2 

3a0 


(42) 


It  may  be  seen  then  that  (at  least  approximately  at  high  Rey¬ 
nolds  number)  ra  scales  with  rv ,  since  according  to  the 
Kolmogorov  hypotheses  a0  is  a  universal  constant.23  In  fact, 
it  is  known11,24-26  that,  at  moderate  Reynolds  numbers,  a0 
increases  weakly  with  Rx — in  accord  with  the  refined  Kol¬ 
mogorov  hypotheses.  In  discussing  scalings  we  ignore  this 
weak  dependence  and  write  ralr—  Re~"I/2. 


The  subsequent  analysis  confirms  that  these  scalings  are  ap¬ 
propriate,  in  that  each  of  these  scaled  coefficients  is  of  order 
unity.  With  these  definitions,  Eq.  (14)  can  be  written 


Tn  \  1/2  ~  «* 


C-A.+  ^  n.ZL 

^IJ  ~  r  ~  -  U l]  ..  f 


dt  _  dWj 

—  +  Bij  - — . 

7 a  dra 


(47) 


Clearly  Ta  is  the  characteristic  time  scale  of  the  process:  the 
mean  of  the  term  in  C,  and  the  variance  of  the  term  in  B  is 
each  of  order  dtfra.  However,  the  term  in  D  is  smaller  by 
the  factor  of  {Tafr)m~ Re“l/4. 

If  the  mean  velocity  gradients  are  constant,  then  the 
equations  for  u*(r)  and  a°(0  [Eqs.  (25)  and  (26)]  can  be 
re-expressed  as  equations  for  u *(r)  and  a *(/)  [Eqs.  (31)  and 
(32)].  The  scaled  forms  of  these  equations  are 


r  duf  /  r\il2af 

^ir=[va )  V 


(48) 


and 
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daf 


~  af  I  ra\  -  af 

c  -  -L-  -j-  I  — —  1  if  J 

ija'  \TSJK^a' 


+ 


1/2 


Dr, 


u * 
J 


dWj 

X 


(49) 


For  the  case  considered  t!ts  is  of  order  unity^so  that  (com¬ 
pared  to  the  leading-order  terms)  the  terms  in  K  and  D  are  of 
order  tJt~  Re“I/2  and  (ra/r)1/2~Re~1/4,  respectively. 


C.  Isotropic  turbulence 

We  consider  in  this  section  the  simplest  case  of  homo¬ 
geneous  isotropic  turbulence  made  statistically  stationary  by 
artificial  forcing.  We  do  so  to  relate  the  general  model  pro¬ 
posed  here  to  Sawford's,12  and  to  provide  a  characterization 
of  the  model's  behavior  in  this  simple  setting.  This  provides 
a  useful  reference  for  the  results  obtained  below  for  the  gen¬ 
eral  case. 

For  isotropic  turbulence  without  mean  velocity  gradi¬ 
ents,  there  is  no  rapid  pressure,  and  a*(t)=du*/dt  is  the 
model  for  the  fluid-particle  acceleration.  The  coefficients  in 
the  model,  Eq.  (32),  are  inevitably  isotropic  (8^=83^, 
Bij  —  SSij ,  etc.),  and  so  the  three  components  of  a *(f)  are 
statistically  identical  and  independent.  Writing  a*(t)  for  one 
component  of  acceleration  [e.g.,  a*(t)^a* (f)],  and  with 
w*(f)  being  the  corresponding  component  of  velocity,  the 
model  for  isotropic  turbulence  is 

da*=-lCa*  +  Du*]dt  +  BdW.  (50) 

This  is  identical  to  Sawford’s  model,12  but  with  the  coeffi¬ 
cients  expresses  differently. 

An  analysis  of  Eq.  (50)  (see  Refs.  12  and  13  and  Sec. 
IV  B)  shows  that  the  acceleration  variance  is 


,a*^EL=a'*t 

'  2 C  °  2C 


(51) 


the  velocity  variance  is 


/  *2v  B  ,2  B 

{u'V  )~  -  •  =  u  — 


2  CD 


2  CD 


(52) 


and  that  the  Lagrangian  velocity  integral  time  scale  is 


where  p(s)  is  the  Lagrangian  velocity  autocorrelation  func¬ 
tion  defined  by  Eq.  (2). 

There  is  a  one-to-one  correspondence  between  the  three 
model  coefficients  8 ,  C,  and  Z),  and  the  three  primary  statis¬ 
tics  a'2,  w'2,  and  TL.  Equations  (5 1)— (53)  are  readily  in¬ 
verted  to  yield  for  the  scaled  coefficients 

C=— ,  B2^—,  0=1.  (54) 

T  T 

The  velocity  autocorrelation  function  p(s)  obtained 
from  the  model  is  most  conveniently  and  naturally  written  in 
terms  of  two  different  (but  related)  time  scales,  and 


r0  (Ta>>r0).  These  are  the  inverses  of  the  two  eigenvalues 
of  the  system,  which  are  given  by  the  solution  to  the 
quadratic  equation 

\2-CX  +  Z>  =  0.  (55) 


The  solutions  are  given  in  terms  of  T L  and  ra  by 


1  + 


4rar\1'2' 

~rTl  . 


and 

1 

X2  =to  =  2Tl 

Conversely  we  have 


Ti-Too+Tq 


and 


T„t0 

T 


(56) 


(57) 

(58) 


(59) 


It  may  be  observed  that  as  ra/TL  tends  to  zero,  and  r0 
tend  to  Tl  and  Ta(riTL),  respectively.  The  coefficients  8 ,  C, 
and  D  given  by  Eq.  (54)  can  be  re-expressed  in  terms  of  T & 
and  r0,  which  is  the  form  originally  given  by  Sawford.12 

The  velocity  autocorrelation  function  given  by  the  model 
is 


p(s)  = 


e 


(60) 


which  is  a  linear  combination  of  two  decaying  exponentials, 
with  time  scales  r0  and  T . 

To  conclude,  based  on  this  examination  of  the  model  in 
isotropic  turbulence,  we  summarize  some  important  observa¬ 
tions  which  are  mirrored  in  the  analysis  of  the  general  model 
presented  below. 

(1)  The  three  model  coefficients  8 ,  C,  and  D  are  uniquely 
related  to  the  three  primary  statistics,  a'2,  w'2,  and  TL. 

(2)  The  autocorrelation  function  p(s)  is  a  linear  combina¬ 
tion  of  decaying  exponentials,  the  time  scales  of  which 
are  the  inverses  of  the  eigenvalues  of  the  system. 

(3)  The  predictions  of  the  model  are  in  excellent  agreement 
with  Lagrangian  statistics  obtained  from  DNS  (see  Refs. 
12  and  13). 

(4)  Given  the  primary  statistics,  a  separate  acceleration  time 
scale  cannot  be  imposed  on  the  model:  instead  the  accel¬ 
eration  time  scale  ra  is  given  by  Eq.  (39). 

(5)  The  simplest  scaling  arguments  show  that  TL  scales  with 
r,  and  that  Ta  scales  with  r^,  so  that  the  scaled  coeffi¬ 
cients  B,  C,  and  D  are  of  order  unity. 


D.  Gaussianity 

For  homogeneous  turbulence,  the  model  takes  the  form 
of  a  set  of  SDEs,  Eqs.  (31)  and  (32),  in  which  the  drift 
coefficients  (-Cuaf  and  -D^uf)  are  linear  in  the  depen¬ 
dent  variables,  while  the  diffusion  coefficient  Btj  is  indepen¬ 
dent  of  a*  and  u*.  Such  linear  stochastic  differential  equa- 
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tions  are  known27  to  yield  Gaussian  processes.  Thus, 
according  to  the  model,  the  processes  a *(0  and  u *(7)  are 
jointly  Gaussian. 

For  homogeneous  turbulent  shear  flow,  the  experiments 
of  Tavoularis  and  Corrsin28  clearly  show  that  the  one-point 
one-time  joint  PDF  of  velocity  is  joint  normal.  Hence  the 
model  is  correct  in  predicting  that  the  one-time  PDF  u*(f)  is 
joint  normal.  However,  it  is  known  from  DNS11  and 
experiments26  29  that  both  acceleration  and  two-time  velocity 
statistics  depart  from  Gaussianity,  an  effect  which  is  not  rep¬ 
resented  by  the  model.  It  is  possible  to  represent  these  effects 
in  stochastic  models  by  making  the  model  coefficients  them¬ 
selves  stochastic  processes.30,31  In  particular.  Beck31  shows 
that  the  experimental  acceleration  distribution  can  be  accu¬ 
rately  represented  by  a  stochastic  model  with  gamma- 
distributed  coefficients.  Here,  however,  we  retain  constant 
coefficients  and  do  not  attempt  to  represent  these  higher- 
order  effects. 

It  is  emphasized  that  the  Gaussianity  of  the  model  is 
confined  to  homogeneous  turbulence.  For  inhomogeneous 
flows,  non-Gaussian  statistics  such  as  the  velocity  triple  cor¬ 
relation  can  be  accurately  calculated  by  linear  stochastic 
models. 


E.  High  Reynolds  number  and  local  isotropy 

We  now  consider  the  limit  of  very  high  Reynolds  num¬ 
ber,  which  is  equivalent  to  the  limit  of  tJt  tending  to  zero. 
In  this  limit,  according  to  the  Kolmogorov  hypotheses,  the 
turbulence  is  locally  isotropic.  As  is  now  shown,  the  stochas¬ 
tic  model  for  acceleration  is  consistent  with  Jocal  isotropy 
provided  that  the  scaled  coefficients  B  and  C  tend  to  the 
following  isotropic  constant  tensors: 

Sy“2(&)-‘*y  and  £,,=(&,)"%.  (61) 

where  C0  is  the  Kolmogorov  constant  associated  with  the 
second-order  Lagrangian  structure  function  [see  Eq.  (69)]. 

In  general,  variations  in  u*(r)  and  a°(t)  occur  on  the 
time  scales  r  and  ra ,  respectively.  For  the  case  considered, 
Ta<r ,  u*(f)  changes  very  slowly  compared  to  a °(r);  and  so 
a °(t)  is  in  a  statistically  quasistationary  state,  the  statistics  of 
which  change  slowly  in  response  to  the  changes  in  u*(r). 
This  state  is  governed  by  Eq.  (47),  with  the  coefficients 
given  by  Eq.  (61),  which  can  be  rewritten 

o  n  dt  a'dWi 

rfao=-(ao_^.[u*(,)])_ - +  -=,  (62) 

4^0  Ta  VfA)ra 

with 


/*<(»*)■- jC0a’  -f  5 y-f  =  -7Co 


1/2 


DijUf 


(63) 


With  ft;  being  considered  as  a  frozen  coefficient,  Eq.  (62)  is 
simply  the  Lange vin  equation;  and  hence  each  component  of 
a°(/)  is  an  independent  Omstein-Uhlenbeck  (OU)  process 
with  conditional  mean  yu,,(u*(/)),  variance  a'2,  and  time 
scale  jC0Ta‘  The  normalized  mean  /x.  /a'  tends  to  zero  as 
(ra/r)  tends  to  zero  [see  Eq.  (63)],  and  hence  a°(/)  tends  to 
a  locally  isotropic  process. 


We  now  examine  the  model  equation  for  velocity  in  the 
high  Reynolds  number  limit.  For  a  general  inhomogeneous 
flow,  the  model  for  u*(/)  [Eqs.  (13)  and  (19)]  is 


duf 

~dT 


(r) 


+g); 


«/+<*<(/>. 


(64) 


As  tJt  tends  to  zero,  a °(f)  tends  to  white  noise;  or,  more 
precisely,  for  a  time  interval  St  such  that  both  raISt  and 
St!  t  tend  to  zero,  the  increment  in  velocity 


(65) 


tends  to  a  Gaussian  random  vector  with  mean  fi(u*(t))St , 
Eq.  (63),  and  covariance 

2a,2QC0Ta)SijSt=C0eSijSt.  (66) 

Thus,  in  the  limit,  Eq.  (64)  tends  to  a  diffusion  process  given 
by  the  SDE, 

du*=[~d-^  +  Gi}\u*dt+(C<)e)mdWi,  (67) 

with 

<  ,  3  Djt 

GiJmGV~4^-y-  (68) 


It  may  be  recognized  that  Eq.  (67)  is  the  generalized 
Langevin  model  (GLM7,9);  and  from  this  observation  we 
draw  two  important  conclusions.  First,  it  is  well  known  that 
the  GLM  is  consistent  with  local  isotropy  and  the  Kolmog¬ 
orov  hypotheses  in  yielding  (for  the  second-order  Lagrangian 
structure  function) 


([wf(/+^)-wf(r)][w;(/+5)-u;(r)]} 

= C0esSu ,  for  s<r. 


(69) 


Second,  in  the  high  Reynolds  number  limit  being  considered, 
Eq.  (68)  gives  the  GLM  coefficient  G;j  which  corresponds  to 
the  acceleration  model  coefficients  G-yr)  and  Di} . 

For  forced,  statistically  stationary  homogeneous  isotro¬ 
pic  turbulence,  the  GLM  coefficient  G(J  is  constrained  to  be 
-jC0Sij/r.9  Correspondingly,  Eq.  (68)  yields  D/;  =  <S,7, 
consistent  with  Sawford's  model,  Eq.  (54).  In  general,  if  the 
GLM  coefficient  is  decomposed  into  slow  and  rapid  con¬ 
tributions,  i.e., 


G-^G^  +  G^ 

lJ  ij  u  ’ 

then  Eq.  (68)  yields 


(70) 


,  ,  3  Du 

rz(s)— _ r _ — 

Gij  ~  4Co  r  • 


(71) 


The  simplest  specification  of  G$j}  (for  unforced  turbulence) 
is 


(72) 


Phys.  Fluids,  Vol.  14,  No.  7,  July  2002 


A  stochastic  Lag  rang  ian  model  for  acceleration  2367 


for  which  the  corresponding  value  of  is 


D 


U 


Sir 


(73) 


In  summary,  with  the  coefficients  B  and  C  specified  by 
Eq.  (61),  the  model  is  consistent  with  the  Kolmogorov  hy¬ 
potheses.  At  very  high  Reynolds  number  (corresponding  to 
tJt  tending  to  zero),  the  acceleration  statistics  are  locally 
isotropic,  and  the  model  tends  to  the  generalized  Langevin 
model  (GLM)  for  velocity,  Eq.  (67).  There  is  then  a  one-to- 
one  correspondence  between  the  remaining  acceleration 
model  coefficient  D  and  the  GLM  coefficient  Gi} ,  Eq.  (68) 
and  Eq.  (71). 


IV.  HOMOGENEOUS  TURBULENT  SHEAR  FLOW 


U*(<) 

»(r)=— ,  - 


„  dt  .  A  d\V(t)  ,  v 

dt=— ,  dW(t)=-ur-,  (11) 

T  T 

these  stochastic  model  equations  transform  to 

dui(t)  =  ai(i)di,  (78) 

dai{i)--iCijaj{i)  +  bijUj(t)‘\dtJt-BijdWj{i),  (79) 

where  the  transformed  (nondimensional)  coefficients  are 


C tj  ~rCij=— C//+  rKi} , 


ij  ■  ’ 


D;;  T^D ij  Djj , 


(80) 


In  this  section  we  examine  the  stochastic  model  for  ac¬ 
celeration  applied  to  homogeneous  turbulent  shear  flow  for 
which  there  are  Lagrangian  data  from  DNS.16’17  The  analysis 
(performed  in  Sec.  IV  B)  depends  on  the  processes  consid¬ 
ered  being  statistically  stationary.  We  therefore  define  (in 
Sec.  IV  A)  a  scaled  time  i,  a  scaled  velocity  u (?),  and  the 
acceleration  a(t)^du(i)ldt  such  that  u(?)  and  a(?)  are  sta¬ 
tistically  stationary  processes — at  least  to  a  reasonable  ap¬ 
proximation.  Results  from  the  analysis  are  compared  to  the 
DNS  data  in  Sec.  IV  C. 

A.  Scaling  for  statistical  stationarity 


1 .  Forced  homogeneous  turbulent  shear  flow 

We  consider  first  the  case  of  forced  homogeneous  turbu¬ 
lent  shear  flow  corresponding  to  the  DNS  of  Schumacher. 14 
This  case  is  relatively  simple  because  the  flow  is  statistically 
stationary.  The  imposed  shear  rate  S  is  constant,  as  are  the 
turbulent  kinetic  energy  k  and  its  dissipation  rate  €.  The  non- 
dimensional  time  t  is  defined  by 


t^t 


€ 

k 


t 

, 

r 


(74) 


and  u (t)  is  defined  as  the  model  for  the  fluctuating  compo¬ 
nent  of  velocity  following  the  fluid  particle,  u*(0,  normal¬ 
ized  by  w ' : 


U  (t)m 


U*(0 

uf 


(75) 


With  these  definitions,  the  velocity  covariance  (ufij)  is  of 
order  unity,  and  so  also  are  the  integral  time  scales  of  u (t) 
(in  scaled  time).  In  fact,  because  of  the  equality  of  one-point, 
one-time  Eulerian  and  Lagrangian  statistics  in  homogeneous 
turbulence,  we  have  the  normalization  condition  following 
from  Eq.  (75): 

(«,(?)«,(?))  =  3.  (76) 

Since  the  velocity  gradients  are  constant,  the  general  stochas¬ 
tic  model  for  u*(f)  and  a*(/)  is  given  by  Eqs.  (31)  and  (32). 
With  the  transformations 


and 

r372  r 

For  a  given  orientation  of  the  shear,  i.e.,  d(  U-)ldxj 
-SSnSJ2,  the  coefficients  B,  C,  and  D  are  constant  and 
depend  only  on  the  Reynolds  number. 


2 .  Unforced  homogeneous  turbulent  shear  flow 

The  DNS  of  Sawford  and  Yeung16  are  consistent  with 
the  supposition  that  (after  an  initial  transient)  the  energy- 
containing  motions  in  (unforced)  homogeneous  turbulent 
shear  flow  become  (approximately)  self-similar.  The  normal¬ 
ized  Reynolds-stress  tensor  ( UjUj)!k  becomes  constant,  as 
does  the  ratio  of  the  turbulence-to-shear  time  scales,  rfrs 
=  Ski  e ,  and  hence  also  the  ratio  of  production  V  to  dissipa¬ 
tion  6.  (The  values  deduced  from  the  DNS  are  Skf  e= 4.83 
and  7Ye=1.54.)  The  turbulent  kinetic  energy  equation  then 
dictates  that  k  and  €  increase  exponentially  with  time — as  is 
observed. 

As  previously  argued,15  this  picture  suggests  that  the 
definitions  of  t  and  u (i)  by  Eq.  (74)  and  Eq.  (75)  remain 
appropriate,  although  now  the  velocity  scale  u'(t)  used  in 
Eq.  (75)  depends  on  time.  This  time  dependence  is  quantified 
by  the  parameter 


r  da’ 

n=^77T 


(82) 


the  value  of  which  is  n~0.27  in  the  present  case.  (The  value 
is  FI=0  for  the  forced  case,  and  11=  -j  for  decaying  turbu¬ 
lence.)  Given  the  (approximately)  self-similar  state  of  the 
energy-containing  motions,  it  is  reasonable  to  suppose  that 
u(f)  is  (approximately)  statistically  stationary.  But  these 
states  can  only  be  realized  approximately  since  the  Reynolds 
number  increases  with  time.  Hence,  while  we  again  define 
a (i)  as  the  derivative  of  u(?),  this  process  cannot  be  com¬ 
pletely  stationary:  according  to  Kolmogorov  scaling,  the  am¬ 
plitude  of  a  increases  as  R {/2  and  its  time  scale  decreases  as 

A  quantification  of  the  variation  of  Rk  in  homogeneous 
turbulent  shear  flow  shows  that  the  departure  from  stationar- 
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ity  is  not  large.  Based  on  the  exponential  increase  of  k  with 
time  it  can  be  shown  that  Rx  increases  as  /?x~exp(II/),  and 
the  DNS  data  are  consistent  with  this  behavior  (except  at  the 
beginning  and  end  of  the  simulation).  The  normalized  La- 
grangian  velocity  integral  time  scale  is  found  to  be  TJt 
=  0.3. 15  Hence,  over  a  time  interval  of  2TL ,  Rx  increases  by 
a  factor  of  exp(0.27X0.6)^1.18.  Thus,  over  the  relevant 
time  interval,  the  amplitude  of  a(?)  increases  by  approxi¬ 
mately  10%,  while  its  time  scale  decreases  by  about  20%. 

As  in  the  forced  case,  the  model  equations  [Eqs.  (31)  and 
(32)]  for  u*(r)  and  a*(0  can  be  transformed  into  equations 
for  u (?)  and  a(f).  The  transformations  are  those  given  by  Eq. 
(77),  except  that  a(0  is  given  by 


a (i)m 


du(t) 

dt 


dfu*(t)\ 

Tdt[u'(t)J 


ra*(0 


— n  u(i). 


(83) 


The  transformed  model  equations  are  again  Eqs.  (78)  and 
(79),  with  B  given  by  Eq.  (81),  but  with  C  and  D  given  by 

Cij~rCij+ 2HSijy 

Du=T2DiJ+TnctJ+n2slj.  (84) 

It  may  be  noted  that  Eq.  (84)  for  C/;-  and  also  applies  to 
the  forced  case,  since  in  that  case  II  is  zero. 

To  conclude,  the  stochastic  model  Eq.  (78)  and  Eq.  (79) 
is  analyzed  in  the  next  section,  with  the  assumptions  that 
u (?)  and  a(?)  are  statistically  stationary.  For  homogeneous 
turbulent  shear  flow,  the  departures  from  stationarity  are  suf¬ 
ficiently  small  that  the  results  of  the  analysis  can  usefully  be 
compared  to  the  DNS  data  of  Sawford  and  Yeung. 17  This  is 
done  in  Sec.  IV  C. 


B.  Analysis  of  the  stochastic  model 

In  this  section  we  analyze  the  model  in  application  to 
homogeneous  turbulent  shear  flow.  The  analysis  is  somewhat 
involved:  for  the  reader  wishing  to  avoid  the  details,  the 
principal  results  are  summarized  in  Sec.  IV  B  6. 


1.  Model  equations 

When  written  for  the  scaled  variables  u(?)  and  a (f)  in 
homogeneous  turbulent  shear  flow,  the  model  equations  are 
Eqs.  (78)  and  (79),  and  the  coefficients  are  given  by  Eqs. 
(81)  and  (84). 

It  is  convenient  to  use  vector-matrix  notation,  and  hence 
we  write  the  model  equations  as 

da(  ?)=-[  Ca(  i)  +  Dfi(  t)  ]di+  B  d W(  i) ,  (85) 

du(i)~a(i)dt,  (86) 

where  the  coefficients  B,  C,  and  D  are  3X3  matrices.  Fur¬ 
thermore,  it  is  convenient  to  combine  a(/)  and  u(/)  into  a 
six-vector 


so  that  the  model  can  be  written  as  the  single  SDE 
dz(i)^-Fz(t)dt  +  Ed\V(t). 


(87) 

(88) 


Here  W(  t )  is  a  six-vector- valued  Wiener  process,  and  the 
6X6  matrix  coefficients  E  and  F  are 


0 

0 


(89) 


and 


C  D 

F=[_,  „]■  m 

where  I  is  the  3X3  identity  matrix. 

It  is  known  from  the  theory  of  diffusion  pro¬ 
cesses927,32’33  that  the  diffusion  coefficient  (e.g.,  B)  affects 
the  process  only  through  the  symmetric  positive-semi- 
definite  from  BBr,  where  ‘T”  denotes  the  transpose.  Hence, 
without  loss  of  generality,  B  and  therefore  E  can  themselves 
be  taken  to  be  symmetric  positive  semidefinite. 

It  is  assumed  that  the  eigenvalues  of  the  drift  matrix  F 
have  positive  real  parts,  which  is  a  sufficient  condition  for 
Eq.  (88)  to  yield  a  statistically  stationary  solution.33 


2.  Autocovariance 


Since  z (/)  is  a  Gaussian  process,  its  statistics  are  com¬ 
pletely  described  by  its  autocovariance,  which  we  define  by 

R(5)=(z(r+s)z(f)7).  (91) 


It  should  be  noted  that  this  is  the  transpose  of  the  conven¬ 
tional  definition  in  that  the  time  increment  s  appears  in  the 
first  variable.  The  present  definition  yields  simpler  equations 
in  the  subsequent  analysis. 

The  autocovariance  of  z (?)  can  be  decomposed  into  the 
autocovariances  of  a(r)  and  u(f): 


R aa(s)  R au(s) 
R  ua(s)  R uu(s) 


(92) 


where 


R “(j)-<u(/+j)a(f)7>,  (93) 

and  Raa(s),  Rau(j),  and  Rww(j)  are  similarly  defined. 

In  view  of  statistical  stationarity,  the  autocovariances  are 
independent  of  time  i  (as  implied  by  the  notation),  and  they 
possess  the  following  symmetries: 

R(s)  =  R(-j)r,  R‘'a(j)  =  R',a(-s)7', 

R"u(s)  =  R"“(- j)r,  (94) 

R"a(i)  =  'Ra"(-*)r==-Ra"(.y).  (95) 

Stemming  from  the  definition  a-duldt,  properties  of 
derivatives  of  the  autocovariances  are 


^Ruu(s)  =  Rau(s),  ^Rua(s)  =  Raa(s), 
d 

—  Rau(s)=-Raa(s). 
ds 


(96) 


and  hence 

<72R',u(s) 

d? 


Rao(j). 


(97) 
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Thus,  all  autocovariances  [including  R(5)]  can  be  deter¬ 
mined  from  RwwO). 

The  covariances  are  denoted  by 


Q=R(0)  = 


Qaa 

Qua 


Qau 

QUU 


(98) 


The  covariance  matrices  Q,  Qaa ,  and  QMM  are  symmetric 
positive  definite;  while  the  off-diagonal  matrices  are  [in  view 
of  Eq.  (95)]  antisymmetric  and  the  transposes  of  each  other, 

QM*  =  ~(Qua)T=(Qau)T.  (99) 

It  may  be  observed  from  Eqs.  (96)  and  (97)  that  all  of  the 
covariances  can  be  obtained  from  RwuO)  and  its  derivatives 
at  the  origin  (5  =  0). 

An  important  quantity  in  the  subsequent  analysis  is  the 
autocorrelation  matrix  which  is  defined  by 


P(s)=R(*)Q-', 


(100) 


and  which  has  the  property 
P(0)  =  I. 


(101) 


3.  Integral  time  scales 

The  matrix  T  of  integral  time  scales,  which  also  plays  a 
central  role  in  the  analysis,  is  defined  by 

T-  f°°P (s)ds.  (102) 

Jo 

This  matrix  has  a  special  structure,  now  revealed,  which 
stems  from  the  fact  that  a  (?)  is  the  derivative  of  u(?).  We 


define  the  6X6  matrix  M  by 

M=  I  R(5)^5, 

Jo 

003) 

which  is  related  to  T  by 

T=MQ“'  or  M=TQ, 

(104) 

T= 


0 


-I 

•J>UU  * 


(108) 


since  the  first  row  of  the  product  TQ  yields  the  first  row  of 
M  [given  by  Eq.  (107)]  in  accordance  with  Eq.  (104). 

By  analogy  to  Eq.  (104),  we  define  the  velocity  integral 
time  scale  tensor  by 

TUU^MUU(QUU)~\  (109) 


which  is  just  (the  transpose  of)  the  time  scale  tensor  that 
arises  in  the  analysis  of  the  stochastic  model  for  velocity.15 
And  we  define  the  (scalar)  Lagrangian  velocity  integral  time 
scale  by  . 


IryiUU 

3  Lii  • 


(110) 


We  see  below  that  the  autocovariance  R(5) — and  there¬ 
fore  all  other  statistics — are  determined  by  the  covariance 
matrix  Q  and  the  time  scale  matrix  T  (as  previously 
shown15).  Because  of  the  special  structure  of  the  model,  the 
information  content  in  Q  and  T  is  less  than  it  appears  at  first 
sight.  Specifically,  the  symmetric  and  nonsymmetric  6X6 
matrices  Q  and  T  can  be  constructed  from  the  3X3  matrices 
Qaa ,  Qu“,  Qw*,  and  Tuu — which  have  an  information  con¬ 
tent  equivalent  to  three  symmetric  and  two  antisymmetric 
3X3  matrices.  It  is  marvelous — although  most  likely 
inevitable — that  this  is  precisely  the  information  content  in 
the  model  coefficients  B,  C,  and  D. 


4.  Solution  for  the  autocorrelation  matrix 

It  is  readily  deduced  from  the  model  equation,  Eq.  (88), 
that  the  autocovariance  satisfies  the  ODE, 

“■  R( 5 )  =  FR( s )  for  55=0.  (Ill) 

ds 

By  post-multiplying  both  sides  of  this  equation  by  Q“l,  we 
find  that  P(5)  [defined  by  Eq.  (100)]  satisfies  the  same  equa¬ 
tion. 


and  which  is  partitioned  as 

[  mt  Mau ' 

[mt  muu  ’ 

For  MT  we  obtain 


M= 


M««=  I  J&.aa(s)ds  —  I  (a (i+s)Z(i)T)ds 
Jo  Jo 

=  |j  a(f+5)J5a(f)r 

=:<[«(00)  —  u(/)]a(f)r> 

=  -<G(f)S(ry>=-<r 


(105) 


(106) 


A  similar  treatment  can  be  applied  to  and  to  show 
that  M  is  given  by 


M= 


'  ~Qua  —Quu 

QUU  MUU  ■ 

It  then  follows  that  T  is  of  the  form 


(107) 


d 

—  P(5)  =  -FP(5)  for  5^0,  (112) 

ds 

with  the  simple  initial  condition  P(0)=I.  The  solution  to  this 
equation  (satisfying  the  initial  condition)  is33 

^  ("l)” 

P(s)  =  exp(— Fs)  =  2  - r-FV  for  s^O,  (113) 

n=0  til 

as  may  be  verified  by  differentiating  with  respect  to  5.  It  has 
been  assumed  that  the  eigenvalues  of  F  have  positive  real 
part,  which  is  a  sufficient  condition  for  exp(-F5)  to  con¬ 
verge  to  zero  as  s  tends  to  infinity. 

The  matrix  F  deduced  from  the  DNS  (in  Sec.  IV  C)  has 
the  simplest  structure — real  positive  eigenvalues,  and  lin¬ 
early  independent  eigenvectors.  In  that  case  F  can  be  decom¬ 
posed  as 

F=  VAV-1,  (114) 

where  the  columns  of  the  6X6  matrix  V  are  the  eigenvectors 
of  F,  and  A  is  the  6X6  diagonal  matrix  of  eigenvalues.  The 
solution  for  P(5),  Eq.  (113),  can  then  be  re-expressed  as 
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P(5')-Vexp(- As)\  1  for  0,  (115) 

showing  that  P(s)  is  a  linear  combination  of  six  decaying 
exponentials,  the  time  scales  of  which  are  the  inverses  of  the 
eigenvalues. 

For  the  general  case,  the  time  scale  matrix  T  [Eq.  (102)] 
is  obtained  as  the  definite  integral  of  the  solution,  Eq.  (113). 
The  indefinite  integral  is 

J  P(s)Js  =  -F~lexp(-Fs),  (116) 

from  which  we  obtain 


T=  f  P(*)<fa  =  F_1. 

Jo 


(117) 


The  6X6  drift  matrix  F  is  defined  in  terms  of  the  3X3  drift 

A  A 

matrices  C  and  D  in  the  stochastic  model  for  acceleration  by 
Eq.  (90).  Given  this  structure  of  F,  it  is  readily  deduced 
(from  the  equation  FF-1  =1)  that  its  inverse  is 


F“l  = 


0 

D1 


-I 

~  i  *  , 

D~  C 


(118) 


which,  according  to  Eq.  (117),  equals  T.  The  first  row  of  F~ 1 
indeed  matches  that  of  T  [Eq.  (108)],  while  equating  the 
elements  of  the  second  rows  yields 

T^-D1,  Twu  =  D-iC,  (119) 

or  conversely 

T>=(Tua)-\  c=(Twa)”1TttU.  (120) 


(The  assumptions  made  about  F  are  sufficient  to  ensure 
that  D  is  nonsingular.) 

The  important  conclusions  are  that  there  is  a  one-to-one 
correspondence  between  the  drift  coefficients  C  and  D  and 
the  time  scale  matrices  Twu  and  Tua  and  that  the  autocorre¬ 
lation  matrix  P(s)  is  explicitly  determined  by  the  drift  coef¬ 
ficients  C  and  D  [through  Eq.  (90),  Eq.  (114),  and  Eq.  (115)]. 
The  autocovariances  are  given  by 


R(s)  =  P(s)Q=exp(  “Fs)Q  for  s^0,  (121) 

where  F  is  given  in  terms  of  C  and  D  by  Eq.  (90). 


5.  Solution  for  the  covariance  matrix 

The  solution  is  completed  by  determining  the  covariance 
matrix  Q.  An  evolution  equation  for  the  covariance  is  readily 
derived  from  the  model  equation  [Eq.  (88)],  and  then  the 
condition  that  Q  is  independent  of  time  yields 

EEr=  E2  =  FQ+  (FQ)r.  (122) 

Thus  E2  is  twice  the  symmetric  part  of  FQV 

From  the  definition  of  E  in  terms  of  B  [Eq.  (89)]  we 
have 

\(E2)aa  (E2)au 

E  _[(E2)ua  (E2)"" 

and  from  the  definitions  of  F  [Eq.  (90)]  and  Q  [Eq.  (92)]  we 
have 


B2  0 

0  0 


(123) 


FQ= 


CQaa  +  I)Qua 

—  Qaa 


CQau+T)Quu 

-Qau 


(124) 


Equation  (122)  can  be  used  to  relate  the  blocks  of  E2  to  FQ, 
and  evidently  [from  Eq.  (123)]  only  the  upper  left-hand 
block  is  nonzero. 

For  the  lower  right-hand  block  we  have,  correctly. 


(E2)um—  —  Qau  —  QauT—Q, 


(125) 


in  view  of  the  antisymmetry  of  Qau  [Eq.  (99)];  and  in  the 
Appendix  it  is  shown  that  the  off-diagonal  blocks  are  also 
zero.  Thus,  the  only  nonzero  block  of  E2  given  by  Eq.  (122) 
is 


( E2y aa  =  B2 - ( CQaa  +  DQua ) +  ( CQaa  +  DQua )T.  (1 26) 


6.  Conclusions 

The  major  conclusion  now  drawn  from  the  analysis  is 
that  there  is  a  one-to-one  correspondence  between  the  model 
coefficients  (B,  C,  and  D)  and  the  primary  statistics  (Q  and 
T).  These  primary  statistics  are  known  in  terms  of  the  veloc¬ 
ity  and  acceleration  covariances  Quu ,  Qua ,  and  Qaa  and  the 
velocity  integral  time  scale  tensor  T““,  Eq.  (109). 

Given  Q  and  T,  the  coefficients  C  and  D  are  given  by 
Eq.  (120),  and  then  B2  is  determined  by  Eq.  (126).  ^ 

Conversely,  given  the  coefficients  (B,  C,  and  D),  T  is 
determined  by  Eq.  (119);  and  the  covariances  are  determined 
by  Eq.  (126)  together  with  the  equation 

Qaa  =  CQau  +  I)Quu.  (127) 

This  equation  is  derived  in  the  Appendix,  where  the  solution 
of  Eq.  (126)  and  Eq.  (127)  for  Q  is  also  discussed.  Together 
these  equations  yield  a  linear  system  which  determines  Q, 
but  unfortunately  an  explicit  solution  is  not  evident. 

Once  both  the  model  coefficients  and  primary  statistics 
are  known,  then  the  autocovariance  given  by  the  model 
R(5,)  =  P(5)Q  can  be  determined  from  Eq.  (115).  These  au¬ 
tocovariances  are  linear  combinations  of  the  six  decaying 
exponentials,  exp(“X^),  where  {\{  Ae}  are  the  ei¬ 

genvalues  of  the  coefficient  matrix  F,  Eq.  (90). 

The  analysis  is  complete,  since  the  autocovariances  Q(s) 
fully  characterize  the  Gaussian  model  processes  a  (?)  and 
u(f). 

For  Sawford's  model  for  isotropic  turbulence,  the 
eigenvalues  of  F  are  t/T «>  and  r/r0  (each  with  multiplicity 
3),  corresponding  to  time  scales  T and  r0  [Eq.  (56)  and 
Eq.  (57)]  which  scale  with  the  integral  time  scale  and 
Kolmogorov  time  scale,  respectively.  And  the  time  scale 
matrices  are 


'UU  _  rjfUU  — 


tq\  ^  Tl  ^ 
r  /  r 


(128) 


and 


'ua 


_T^tq 


1= 


T 


(129) 
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FIG.  1.  Velocity  autocovariances 
R ““(s)  against  normalized  time.  Sym¬ 
bols,  DNS  of  Sawford  and  Yeung 
(Ref.  16);  lines,  from  the  acceleration 
model. 


C.  Comparison  to  DNS  data 

In  this  section,  for  the  DNS  of  homogeneous  turbulent 
shear  flow,16  the  stochastic  model  coefficients  B,  C,  and  D 
are  deduced  from  the  data;  and  then  the  velocity  and  accel¬ 
eration  autocovariances  predicted  by  the  model  are  compared 
to  those  from  the  DNS. 

All  the  DNS  information  is  extracted  from  the  time  se¬ 
ries  of  the  normalized  velocity  autocovariance  Rww(j).  The 
remaining  autocovariances  [RMa(5),  Ra“(^),  and  Raa(s)]  are 
obtained  from  Eq.  (96)  and  Eq.  (97)  by  numerical  differen¬ 
tiation  of  Ruu(s),  and  then  the  covariances  are  obtained  as 
Q==R(0).  Clearly  this  differentiation  amplifies  the  statistical 
noise  in  the  data,  as  is  particularly  evident  in  Raa(s)  (see 
Fig.  3  below).  [In  future  DNS  designed  for  this  purpose,  it 
would  be  preferable  to  form  all  covariances  directly  from 
a  (?)  and  u(?).] 

The  velocity  covariance  integrals  [Eqs.  (103)  and 
(105)]  are  formed  from  the  time  series  of  Rw“(s)  by  numeri¬ 
cal  quadrature,  and  then  the  matrix  of  integral  time  scales  T 
is  obtained  from  Eq.  (104). 

Based  on  the  DNS  values  of  the  covariances  and  the 
velocity  integral  time  scales,  the  values  of  the  model  co¬ 
efficient  B,  C,  and  D  are  deduced  which  lead  to  the  model's 
matching  these  statistics.  The  values  of  C  and  D  are  obtained 
from  Eq.  (120),  and  then  B2  from  Eq.  (126).  (The  values  thus 
obtained  are  reported  below.)  The  autocovariances  predicted 
by  the  model  are  then  deduced  from  Eqs.  (114),  (115),  and 
(100). 

Figures  1-3  show  a  comparison  of  the  autocovariances 
from  the  DNS  (symbols)  and  from  the  model  (solid  lines). 
Clearly  the  agreement  is  excellent,  especially  for  the  velocity 
autocovariances  (Fig.  1).  It  should  be  noted  that  an  accelera¬ 
tion  time  scale  is  not  an  input  to  the  model,  and  so  the 


matching  of  the  location  and  magnitudes  of  the  peaks  of 
R“a(s)  and  Raa(s)  in  Figs.  2  and  3  is  not  inevitable. 

Figure  4  compares  the  velocity  autocovariances  from  the 
DNS,  from  the  present  acceleration  model,  and  from  the  sto¬ 
chastic  model  for  velocity15  (dashed  lines).  Only  the  early 
times  are  shown  where  the  differences  between  the  two  mod¬ 
els  are  most  evident.  As  may  be  seen,  the  acceleration  model 
provides  a  much  more  accurate  representation  of  the  curva¬ 
ture  of  the  autocovariances  at  small  times.  As  the  Reynolds 
number  increases,  the  differences  between  the  two  models 
decreases,  and  is  confined  to  smaller  times. 

The  values  of  the  coefficients  B,  C,  and  D  deduced  from 
the  DNS  are  reported  in  scaled  form  as  B  and  D  [defined 
by  Eq.  (46)]  and  ra C.  Their  values  are 


S2  = 


0.70  -0.15  0 

-0.15  0.48  0 

0  0  0.64. 


c= 


0.33 

0.07 

0 


-0.03  0 

0.29  0 

0  0.27 


D= 


0.70  0.15 
0.33  1.41 

0  0 


0 

0 

1.11 


(130) 


(131) 


(132) 


Given  that  d{Ul)/dx2  is  the  only  nonzero  velocity  gradient, 
the  symmetries  in  the  problem  dictate  that  the  off-diagonal 
components  in  the  third  rows  and  columns  of  these  matrices 
are  zero — as  is  observed.  The  magnitudes  of  all  three  coef- 
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FIG.  2.  Velocity-acceleration  auto¬ 
covariances  R  1j(s)  against  normalized 
time.  Symbols,  DNS  of  Sawford  and 
Yeung  (Ref.  16);  lines,  from  the  ac¬ 
celeration  model.  [Note  that  R ™(s) 


ficients  are  as  expected  from  Eq.  (54)  given  that  the  velocity 
integral  time  scale  is  Tl/t^ 0.3. 

As  discussed  in  Sec.  Ill  E,  if  local  isotropy  prevailed  at 
high  Reynolds  number,  then  the  acceleration  statistics  would 
be  isotropic  (to  leading  order  in  *).  A  sufficient  condition 
for  the  model  to  yield  such  local  isotropy  is  that  B  and  C  (but 
not  C)  become  isotropic  as  the  Reynolds  number  increases. 


It  is  evident  from  Eq.  (130)  that,  at  the  moderate  Reynolds 
number  of  the  DNS,  B2  exhibits  significant  anisotropy. 

V.  APPLICATION  TO  TURBULENCE  MODELLING 

The  general  model  proposed  here  consists  of  an  ODE 
Eq.  (11)  for  the  fluid-particle  velocity  U*(/),  which  includes 


FIG.  3.  Acceleration  autocovariances 
R °f(s)  against  normalized  time.  Sym¬ 
bols,  DNS  of  Sawford  and  Yeung 
(Ref.  16);  lines,  from  the  acceleration 
model. 
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a  rapid-pressure  model  [Eq.  (15)];  and  an  SDE  [Eq.  (12)]  for 
the  acceleration  variable  A°(/).  A  specific  model  consists  of 
a  specification  of  the  coefficients  appearing  in  these  equa¬ 
tions,  namely,  Ntjk ,  ,  and  . 

Beyond  proposing  the  general  model,  the  objective  here 
is  not  to  suggest  a  specific  model,  but  rather  to  show  that  all 
of  the  coefficients  can  be  deduced  from  DNS  data  on  homo¬ 
geneous  turbulence.  Hence,  future  DNS  studies — at  different 
Reynolds  numbers  and  with  different  imposed  mean  velocity 
gradients — can  be  used  to  guide  the  construction  of  a  specific 
model. 

As  outlined  in  the  following  subsections,  the  accelera¬ 
tion  model  can  be  used  at  different  levels  of  turbulence  mod¬ 
eling.  In  each  case,  the  turbulent  time  scale  r  is  needed, 
which  can  be  obtained  from  the  standard  model  equation  for 
e  or  o)=  t~  1 ,  or  from  particle  models  for  such  quantities.30,34 

A.  Velocity-acceleration-wave-vector  model 

In  addition  to  the  model  equations  for  U*(r)  and  A°(/), 
an  additional  SDE  can  be  solved  for  the  unit  wave  vector 
e*(f)  (Refs.  9  and  21),  so  that  Eq.  (18)  can  be  used  as  the 
rapid-pressure  model.  Such  a  model  has  the  virtue  of  repre¬ 
senting  exactly  the  evolution  of  the  Reynolds  stresses  for 
arbitrary  rapid  distortions  of  homogeneous  turbulence.  Al¬ 
though  it  has  not  been  convincingly  demonstrated,  the  model 
should  also  be  capable  of  providing  a  more  accurate  repre¬ 
sentation  of  the  rapid  pressure  away  from  the  rapid-distortion 
limit. 

B.  Velocity-acceleration  model 

Without  the  wave-vector  model,  the  rapid  pressure  has  to 
be  modeled  in  terms  of  the  particle  velocity  and  Reynolds 
stresses  (among  other  quantities).  The  standard  model  [Eq. 


FIG.  4.  Velocity  autocovariances 
R““0)  at  early  times.  Symbols,  DNS 
data  of  Sawford  and  Yeung  (Ref.  17); 
solid  line,  from  the  acceleration 
model;  dashed  line,  from  the  velocity 
model  (Ref.  15). 


(19)]  is  linear  in  the  velocity,  but  nonlinear  models  [e.g.,  Eq. 
(24)]  can  also  be  considered.  It  has  to  be  acknowledged  that, 
at  this  level  of  closure,  there  is  insufficient  directional  infor¬ 
mation  to  model  accurately  rapid  distortions.  But  such  mod¬ 
els  may  be  adequately  accurate  for  the  moderate  distortions 
that  typically  occur  in  turbulent  shear  flows. 

Compared  to  a  velocity  model  (discussed  in  the  next 
section),  a  velocity-acceleration  model  has  two  advantages. 
First,  in  essence  it  models  the  velocity  as  a  second-order 
system,  Eq.  (29),  rather  than  as  a  first-order  system.  Conse¬ 
quently,  the  rapid  and  slow  responses  of  the  turbulence  to  a 
sudden  change  in  the  mean  velocity  gradients  can  be  mod¬ 
eled  in  a  natural  way.  The  second  advantage  is  that  accelera¬ 
tion  is  modeled  realistically  rather  than  as  white  noise,  and 
thereby  Reynolds  number  effects  can  be  incorporated  in  a 
natural  way. 

An  apparent  disadvantage  is  that,  in  a  numerical  imple¬ 
mentation,  time  steps  At  of  order  ra  (or  equivalently  tv)  are 
needed  to  resolve  the  acceleration  time  series;  whereas  with 
a  velocity  model  the  time  steps  can  be  of  order  TL  (or 
equivalently  rj*.  With  time  steps  of  order  tv,  the  computa¬ 
tional  cost  increases  as  ReI/2.  However,  in  most  applications 
the  details  of  the  short-time  behavior  are  not  required,  and 
temporal  resolution  on  a  time  scale  of  order  r  is  sufficient.  It 
is  fortunate,  therefore,  that  the  model  equations  can  be 
solved  accurately  by  numerical  methods  that  take  time  steps 
At  that  are  large  compared  to  rv  (but  small  compared  to  r). 
This  is  because  the  model  coefficients  vary  on  the  time  scale 
r,  and  the  model  equations  [e.g.,  Eq.  (25)  and  Eq.  (26)]  with 
frozen  coefficients  admit  analytic  solutions.  Consequently,  if 
resolution  on  the  Kolmogorov  time  scale  is  not  required,  the 
velocity-acceleration  model  can  be  implemented  with  a  com¬ 
putational  cost  that  is  independent  of  Reynolds  number. 
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C.  Velocity  model 

Given  a  specific  velocity-acceleration  model,  a  corre¬ 
sponding  velocity  model  can  be  defined — as  now  outlined. 

When  applied  to  (approximately)  self-similar  homoge¬ 
neous  turbulence  (at  a  given  Reynolds  number  and  with 
given  imposed  mean  velocity  gradients),  the  velocity- 
acceleration  model  yields  a  value  of  the  normalized  Rey¬ 
nolds  stress  tensor  ( Quu )  and  of  the  velocity  integral  time 
scale  tensor  (TWM).  The  corresponding  velocity  model  is  de¬ 
fined  to  be  the  linear  SDE  for  velocity  that  yields  these  same 
statistics.  The  drift  and  diffusion  coefficients  in  the  velocity 
model  are  uniquely  determined  by  Quu  and  Twu.15 

[This  procedure  for  determining  the  velocity-model  co¬ 
efficients  is  straightforward  to  implement  numerically;  but  an 
analytical  treatment  is  hampered  by  the  lack  of  an  explicit 
solution  to  Eq.  (122)  for  Q.] 

The  form  of  the  velocity  model  thus  obtained  is  the  same 
as  the  generalized  Langevin  model  (GLM7,9)  but  with  an 
anisotropic  diffusion  coefficient.  The  advantage  of  obtaining 
a  velocity  model  by  this  route  is  that  it  inherits  the  Reynolds- 
number  dependence  (and  other  attributes)  of  the  velocity- 
acceleration  model.  At  very  high  Reynolds  number  the  ac¬ 
celeration  model  tends  to  the  GLM  with  isotropic  diffusion, 
Eq.  (67),  and  with  the  coefficient  G/y  given  by  Eq.  (71). 

D.  Reynolds-stress  model 

Given  a  particle  model  for  velocity,  it  is  straightforward 
to  derive  a  corresponding  Reynolds-stress  equation. 1020,35 
Again,  such  a  model  inherits  from  its  antecedents  a 
Reynolds-number  dependence  and  other  attributes. 


of  acceleration  and  of  multitime  velocity  statistics  is  physi¬ 
cally  incorrect,  and  reflects  the  fact  that  the  model  does  not 
account  for  internal  intermittency. 

For  homogeneous  turbulent  shear  flow,  the  model  coef¬ 
ficients  are  evaluated  from  the  DNS  data  of  Sawford  and 
Yeung.17  The  model  autocovariances  thus  obtained  (Figs. 
1-4)  are  in  excellent  agreement  with  those  from  the  DNS, 
including  the  short-time  (Kolmogorov  scale)  behavior  (see 
Fig.  4). 

Compared  to  a  linear  stochastic  model  for  velocity,  the 
velocity-acceleration  model  has  the  advantage  of  providing  a 
realistic  representation  of  the  behavior  on  the  Kolmogorov 
time  scale;  and,  as  a  consequence,  of  naturally  incorporating 
Reynolds-number  effects.  The  purpose  here  has  not  been  to 
propose  a  specific  model  (i.e.,  a  specification  of  the  model 
coefficients),  but  rather  to  show  that  these  coefficients  can  be 
deduced  from  DNS  of  homogeneous  turbulence,  as  functions 
of  the  Reynolds  number  and  of  the  imposed  mean  velocity 
gradients. 

As  discussed  in  Sec.  V,  the  velocity-acceleration  model 
can  be  used  as  a  basis  for  generating  a  range  of  (Reynolds- 
number  dependent)  PDF  and  Reynolds-stress  turbulence 
models. 
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VI.  CONCLUSIONS 


APPENDIX:  SOLUTION  FOR  THE  COVARIANCE 
MATRIX 


We  have  considered  a  stochastic  model  for  fluid-particle 
velocity  and  acceleration  in  inhomogeneous  turbulent  flows. 
The  model  consists  of  an  ODE  for  velocity,  Eq.  (11),  and  an 
SDE  for  an  acceleration  variable,  Eq.  (12).  This  structure 
produces  the  correct  qualitative  response  to  rapid  distortions. 
If  the  model  is  supplemented  by  the  wavevector  equation, 
then  the  resulting  model  [Eq.  (15)  and  Eq.  (18)]  is  exact  for 
arbitrary  rapid  distortions  of  homogeneous  turbulence.  Oth¬ 
erwise,  a  standard  linear  model,  Eq.  (19),  for  the  rapid  pres¬ 
sure  can  be  used. 

For  isotropic  turbulence,  the  SDE  for  acceleration  re¬ 
duces  to  Sawford's  model.12  For  very  high  Reynolds 
number  the  model  is  consistent  with  local  isotropy  and 
the  Kolmogorov  hypotheses,  and  tends  to  the  generalized 
Langevin  model  for  velocity.  For  homogeneous  turbulence 
(with  constant  and  uniform  imposed  mean  velocity  gradi¬ 
ents)  a  full  analysis  of  the  model  is  performed.  This  estab¬ 
lishes  the  one-to-one  correspondence  between  the  model  co¬ 
efficient  tensors  B,  C,  and  D  and  the  primary  statistics  of  the 
model,  namely,  the  velocity-acceleration  covariances  and  the 
velocity  integral  time  scale  tensor.  Details  are  given  in  Sec. 
IV  B  6.  For  homogeneous  turbulence,  the  modeled  processes 
(i.e.,  the  velocity  and  acceleration  time  series)  are  Gaussian, 
and  hence  are  completely  characterized  by  their  autocovari¬ 
ance,  which  is  given  explicitly  by  Eq.  (121).  The  Gaussianity 


The  purposes  of  this  appendix  are  to  derive  Eq.  (127);  to 
show  that  the  off-diagonal  blocks  of  E2  given  by  Eq.  (122) 
and  Eq.  (124)  are  zero;  and  to  discuss  the  solution  of  Eq. 
(126)  and  Eq.  (127)  for  the  covariances. 

The  covariances  are  related  to  the  derivatives  of  the  au¬ 
tocovariances  at  the  origin,  Eq.  (96).  From  the  ODE  for  the 
model  autocovariance  [Eq.  (Ill)]  we  obtain 


d 

C  D 

'Qaa 

Qau 

=— FQ=- 

5  =  0 

-I  0 

Qua 

QUU 

-  CQaa  ~  f)Qua  -  CQau  “  T>QUU 

Qaa  quu 

(Al) 

The  bottom  row  of  this  last  matrix  is  consistent  with  the  first 
two  relations  in  Eq.  (96);  while  the  consistency  of  the  upper 
right  block  with  the  third  relation  in  Eq.  (96)  yields 

Qaa=CQau+I)Qu\  (A2) 


which  is  Eq.  (127). 

The  second  derivative  at  the  origin  is 
72 


^R(S) 


=  F2Q. 


5  =  0 


(A3) 
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By  expanding  the  right-hand  side  and  invoking  Eq.  (97)  we 
obtain 


=  Qaa  =  CQau  +  i)Quu, 

s~  0 


(A4) 


which  provides  no  new  information,  but  is  consistent  with 
Eq.  (A2). 

From  Eq.  (122)  and  Eq.  (124),  we  obtain  for  the  upper 
right  block  of  E2, 


CQau  +  DQuu-Qa 


=  Qaa  —  Qa“r=0, 


(A6) 


where  the  second  line  follows  from  Eq.  (A2)  and  the  sym¬ 
metry  of  Qaa . 

We  turn  now  to  the  solution  of  Eq.  (126)  and  Eq.  (127) 
for  Quu ,  Qa*\  and  Qua  given  B,  C,  and  D.  Recall  that 
and  are  symmetric,  while  Qua  is  antisymmetric.  Both 
sides  of  Eq.  (126)  are  identically  symmetric,  whereas  the 
right  hand  side  of  Eq.  (127)  is  not  identically  symmetric. 
Hence,  together,  these  equations  represent  a  linear  system  for 
the  components  of  Qwu,  Qafl,  Qua — with  the  same  number  of 
independent  equations  as  the  number  of  independent  un¬ 
knowns,  i.e.,  15. 

It  is  very  unfortunate  that  there  appears  not  to  be  a 
simple  explicit  solution  for  the  covariances.  It  should  be  pos¬ 
sible,  however,  to  obtain  an  explicit  solution  using  tensor 
representation  theorems.36'37  That  is,  the  covariance  can  be 
'written 


Ns  Ns 

Q"“=X  r(unj s<">,  cr=2  r!")s(n) 

n=  1  n= 1 

Na 

Q“a=2  rinJ A<">, 

n=  1 


(A7) 


where  {S(n>}  is  a  complete  set  of  Ns  linearly  independent 
symmetric  tensor  functions  that  can  be  formed  for  B,  C,  and 
D;  and  similarly  {A^}  is  a  complete  set  of  Na  antisymmet¬ 
ric  tensors.  The  coefficients  {r^},  {r^},  and  can  then 
be  deduced  from  Eq.  (126)  and  Eq.  (127):  they  are  invariants 
of  B,  C,  and  D.  However  such  a  solution  is  unlikely  to  be 
simple  (or  easy  to  obtain). 
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